text:
__precompile__()

module ConjugatePriors

using Statistics
using LinearAlgebra

using PDMats
using Distributions
using StatsFuns
using SpecialFunctions

import Statistics: mean
import LinearAlgebra: Cholesky

import PDMats: PDMat

import StatsFuns:
    logmvgamma,
    log2π

import Distributions:
    BernoulliStats,
    BinomData,
    BinomialStats,
    CategoricalData,
    CategoricalStats,
    DirichletCanon,
    ExponentialStats,
    IncompleteFormulation,
    MultinomialStats,
    MvNormalStats,
    MvNormalKnownCovStats,
    NormalStats,
    NormalKnownSigma,
    NormalKnownMu,
    NormalKnownSigmaStats,
    NormalKnownMuStats,

    add!,
    add_categorical_counts!,
    shape,
    scale,
    rate,
    mode,
    rand,
    pdf,
    logpdf,
    params,
    isApproxSymmmetric,
    hasCholesky

export
    # conjugate prior distributions defined here
    NormalInverseChisq,
    NormalGamma,
    NormalInverseGamma,
    NormalWishart,
    NormalInverseWishart,

    # inteface
    posterior,
    posterior_canon,
    posterior_rand,
    posterior_rand!,
    posterior_randmodel,
    posterior_mode,
    fit_map

include("fallbacks.jl")
include("beta_binom.jl")
include("dirichlet_multi.jl")
include("gamma_exp.jl")

include("normalgamma.jl")
include("normalinversegamma.jl")
include("normalinversechisq.jl")
include("normalwishart.jl")
include("normalinversewishart.jl")
include("normal.jl")
include("mvnormal.jl")

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FSPS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FSPS")
JLLWrappers.@generate_main_file("FSPS", UUID("f5824f2b-013a-5983-a57b-cbb6554a033a"))
end  # module FSPS_jll

ext: jl


text:
module test_01_OneScalarSignal

using ModiaPlot
include("test_base/test_01_OneScalarSignal.jl")

end


text:
#############################
# Interface for distributions
#############################

"""
    rand{T}(rng::AbstractRNG, d::Any)

Return a random element from distribution or space `d`.

If `d` is a state or transition distribution, the sample will be a state; if `d` is an action distribution, the sample will be an action or if `d` is an observation distribution, the sample will be an observation.
"""
function rand end

"""
    pdf(d::Any, x::Any)

Evaluate the probability density of distribution `d` at sample `x`.
"""
function pdf end

"""
    mode(d::Any)

Return the most likely value in a distribution d.
"""
function mode end

"""
    mean(d::Any)

Return the mean of a distribution d.
"""
function mean end

"""
    support(d::Any)

Return the possible values that can be sampled from distribution d. Values with zero probability may be skipped.
"""
function support end

"""
    iterator(d::Any)

DEPRECATED. Return an iterable object (array or custom iterator) that iterates over possible values of distribution `d`. Values with zero probability may be skipped.
"""
@generated function iterator(x::Any)
    @warn("POMDPs.iterator(x) is deprecated. Simply iterate over the space directly or use support(d) for distributions.")
    return :(support(x))
end

"""
    sampletype(T::Type)
    sampletype(d::Any) = sampletype(typeof(d))

Return the type of objects that are sampled from a distribution or space `d` when `rand(rng, d)` is called.

The distribution writer should implement the `sampletype(::Type)` method for the distribution type, then the function can be called for that type or for objects of that type (i.e. the `sampletype(d::Any) = sampletype(typeof(d))` default is provided).
"""
function sampletype end

sampletype(d::Any) = sampletype(typeof(d))
sampletype(t::Type) = throw(MethodError(sampletype, (t,)))

implemented(f::typeof(sampletype), TT::Type{Tuple{T}}) where T<:Type = hasmethod(f, TT) && which(f, TT).module != POMDPs
implemented(f::typeof(sampletype), ::Type{Tuple{T}}) where T = implemented(f, Tuple{Type{T}})


text:
using RBNF
using MLStyle
using PrettyPrint
using DataStructures
struct ReML
end

second((a, b)) = b
second(vec::V) where V <: AbstractArray = vec[2]

escape(a) = @match a begin
    '\\'  => '\\'
    '"'   => '"'
    'a'   => '\a'
    'n'   => '\n'
    'r'   => '\r'
    a     => throw(a)
end

join_token_as_str = xs -> join(x.value for x in xs)

maybe_to_bool(::Some{Nothing}) = false
maybe_to_bool(::Some{T}) where T = true

"""
nonempty sequence
"""
struct GoodSeq{T}
    head :: T
    tail :: Vector{T}
end

join_rule(sep, x) = begin
    :([$x, ([$sep, $x] % second){*}] % x -> GoodSeq(x[1], x[2]))
end

RBNF.typename(::Type{ReML}, name:: Symbol) = Symbol(:R, name)

RBNF.@parser ReML begin
    # define the ignores
    ignore{space}

    # define keywords
    reserved = [true, false]

    @grammar
    # necessary
    Str       :=  value=str

    Bind      := [name=id, '=', value=Exp]
    Let       := [hd=:let, rec=:rec.? % maybe_to_bool,
                        binds=join_rule(',', Bind),
                     :in, body=Exp]

    Fun       := [hd=:fn, args=id{*}, "->", body=Exp]
    Import    := [hd=:import, paths=join_rule(',', id)]
    Match     := [hd=:match, sc=Exp, :with,
                        cases=(['|', Comp, "->", Exp] % ((_, case, _, body), ) -> (case, body)){*},
                    :end.?] # end for nested match
    If        := [hd=:if, cond=Exp, :then,
                        br1=Exp,
                     :else,
                        br2=Exp]
    Num       := [[neg="-"].? % maybe_to_bool, (int=integer) | (float=float)]
    Boolean   := value=("true" | "false")
    NestedExpr = ['(', value=Exp, ')'] => _.value
    Var       := value=id
    Block     := [hd=:begin, stmts=Stmt{*}, :end]
    Atom      =  NestedExpr | Num | Str | Boolean | Var
    Annotate  := [value=Atom, [':', Exp].?]
    Call      := [fn=Annotate, args=Annotate{*}]
    Comp      := [[forall=Forall, '.'].?, value=(Call | Let | Fun | Match | If | Block)]
    Op        := ['`', name=_, '`'] | [is_typeop :: Bool = "->" => true]
    Exp       := [hd=Comp, tl=[Op, Comp]{*}]

    Decl      := [hd=:val, name=id, ':', typ=Exp]
    Define    := [hd=:def, name=id, '=', value=Exp]
    Stmt      = Data | Import | Class | Instance | Decl |  Define | Exp
    Module     := [hd=:module, name=id, :where, stmts=Stmt{*}, :end.?]

    # forall
    Constaints:= [hd=:that, elts=join_rule(',', Call)]
    Forall    := [:forall, fresh=id{*}, (constraints=Constaints).?]

    # advanced
    Data      := [hd=:data, name=id, ":", kind=Exp, :where, stmts=Stmt{*}, :end]
    Class     := [hd=:class, name=id, ids=id{*},  (constrains=Constaints).?, :where,
                     interfaces=Decl{*},
                 :end]
    Instance  := [hd=:instance, name=id, vars=Exp{*}, :where, impls=Stmt{*}, :end]

    @token
    id        := r"\G[A-Za-z_]{1}[A-Za-z0-9_]*"
    str       := @quote ("\"" ,"\\\"", "\"")
    float     := r"\G([0-9]+\.[0-9]*|[0-9]*\.[0.9]+)([eE][-+]?[0-9]+)?"
    integer   := r"\G([1-9]+[0-9]*|0)"
    space     := r"\G\s+"
end

src1 = """
module Poly where
def a = "12\\"3"

class Monad m that Functor m where
    val bind : forall a b . m a -> (a -> b) -> m b
end

val a : Int
def a = "12345"

val f : Int -> Int
def f = fn a -> a `add` 1

val g : Int
def g = match 1 with
        | 1 -> 2
        | _ -> 0
        end

def z = fn x -> if x `equals` 1 then 1 else 2
"""

tokens = RBNF.runlexer(ReML, src1)
ast, ctx = RBNF.runparser(Module, tokens)
pprint(ast)


text:
#This is an odd place to keep these.  Could see about PR to AxisArrays?
axisspacing(A::AxisArray) = map(step, axisvalues(A))
match_axisspacing(B::AbstractArray{T,N}, A::IM) where {T,N, IM<:ImageMeta} = ImageMeta(match_axisspacing(B, data(A)), properties(A))
function match_axisspacing(B::AbstractArray{T,N}, A::AxisArray{T2,N}) where {T,T2,N}
    sp = axisspacing(A)
    nms = axisnames(A)
    axf = ax->begin
        u = unit(sp[ax])
        Axis{nms[ax]}(range(0.0*u, stop=(size(B,ax)-1)*sp[ax], length=size(B,ax)))
    end
    newaxs = map(axf, 1:length(sp))
    return AxisArray(B, newaxs...)
end


text:
using SafeTestsets, Test

@safetestset "Model" begin include("model.jl") end
@safetestset "DictModel" begin include("dict_model.jl") end
@safetestset "JuMP" begin include("jump.jl") end
@safetestset "Functions" begin include("functions.jl") end
@safetestset "Flatten" begin include("flatten.jl") end


text:
@forwardfun bfmi
@forwardfun ess
@forwardfun geweke
@forwardfun mcse
@forwardfun rhat

function convert_arguments(::typeof(bfmi), data, args...; kwargs...)
    dataset = convert_to_dataset(data; group = :sample_stats)
    return tuple(dataset, args...), kwargs
end
function convert_arguments(::typeof(bfmi), data::AbstractArray, args...; kwargs...)
    return tuple(data, args...), kwargs
end

for f in (:ess, :mcse, :rhat)
    @eval begin
        function convert_arguments(::typeof($(f)), data, args...; kwargs...)
            dataset = convert_to_dataset(data; group = :posterior)
            return tuple(dataset, args...), kwargs
        end
        function convert_arguments(::typeof($(f)), data::AbstractArray, args...; kwargs...)
            return tuple(data, args...), kwargs
        end
    end
end


text:
function Get_library_version()
    # There is no way to query at runtime what the length of the buffer should be.
    # https://github.com/mpi-forum/mpi-issues/issues/159
    # 8192 is the maximum value of MPI_MAX_LIBRARY_VERSION_STRING across known
    # implementations.
    buf = Array{UInt8}(undef, 8192)
    buflen = Ref{Cint}()

    libfilename, = split(basename(libmpi),'.')
    @mpicall ccall((:MPI_Get_library_version, libmpi), Cint, (Ptr{UInt8}, Ref{Cint}), buf, buflen)
    resize!(buf, buflen[])
    return String(buf)
end

"""
    MPI_LIBRARY_VERSION_STRING :: String

The full version string provided by the library

# External links
$(_doc_external("MPI_Get_library_version"))
"""
const MPI_LIBRARY_VERSION_STRING = Get_library_version()

"""
    MPIImpl

An enum corresponding to known MPI implementations

- `UnknownMPI`: unable to determine MPI implementation
- `MPICH`: [MPICH](https://www.mpich.org/)
- `OpenMPI`: [Open MPI](https://www.open-mpi.org/)
- `MicrosoftMPI`: [Microsoft MPI](https://docs.microsoft.com/en-us/message-passing-interface/microsoft-mpi)
- `IntelMPI`: [Intel MPI](https://software.intel.com/en-us/mpi-library)
- `SpectrimMPI`: [IBM Spectrum MPI](https://www.ibm.com/us-en/marketplace/spectrum-mpi)
- `MVAPICH`: [MVAPICH](http://mvapich.cse.ohio-state.edu/)

# See also

- [`MPI_LIBRARY`](@ref)
"""
@enum MPIImpl begin
    UnknownMPI
    MPICH
    OpenMPI
    MicrosoftMPI
    IntelMPI
    IBMSpectrumMPI
    MVAPICH
end

"""
    MPIABI

An enum corresponding to known MPI Application Binary Interfaces (ABI)

- `UnknownABI`: unable to determine MPI ABI. MPI.jl will attempt to build a small C
  program to determine the necessary constants and type information.

- `MPICHABI`: Compatible with [MPICH ABI Compatibility
  Initiative](https://www.mpich.org/abi/).

- `OpenMPIABI`: Compatible with [Open MPI](https://www.open-mpi.org/).

- `MicrosftMPIABI`: Compatible with [Microsoft MPI](https://docs.microsoft.com/en-us/message-passing-interface/microsoft-mpi).
"""
@enum MPIABI begin
    UnknownABI
    MPICHABI
    OpenMPIABI
    MicrosoftMPIABI
end

"""
    impl, version, abi = identify_implementation()

Attempt to identify the MPI implementation based on
[`MPI_LIBRARY_VERSION_STRING`](@ref). Returns a triple of values:

- `impl`: a value of type [`MPIImpl`](@ref)
- `version`: a `VersionNumber` of the library, or `nothing` if it cannot be determined.
- `abi`: a value of [`MPIABI`](@ref). This can be overridden by the `JULIA_MPI_ABI` [environment variable](@ref environment_variables).

This function is only intended for internal use. Users should use [`MPI_LIBRARY`](@ref),
[`MPI_LIBRARY_VERSION`](@ref) or [`MPI_LIBRARY_ABI`](@ref).
"""
function identify_implementation()
    impl = UnknownMPI
    version = nothing
    abi = UnknownABI

    if startswith(MPI_LIBRARY_VERSION_STRING, "MPICH")
        impl = MPICH
        # "MPICH Version:\t%s\n" /  "MPICH2 Version:\t%s\n"
        if (m = match(r"^MPICH2? Version:\t(\d+.\d+.\d+\w*)\n", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
            if version >= v"3.1"
                abi = MPICHABI
            end
        end

    elseif startswith(MPI_LIBRARY_VERSION_STRING, "Open MPI")
        # Open MPI / Spectrum MPI
        impl = occursin("IBM Spectrum MPI", MPI_LIBRARY_VERSION_STRING) ? IBMSpectrumMPI : OpenMPI

        if (m = match(r"^Open MPI v(\d+.\d+.\d+\w*)", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
        end
        abi = OpenMPIABI

    elseif startswith(MPI_LIBRARY_VERSION_STRING, "Microsoft MPI")
        impl = MicrosoftMPI
        # "Microsoft MPI %u.%u.%u.%u%S"
        # ignore last 2 (build numbers)
        if (m = match(r"^Microsoft MPI v(\d+.\d+)", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
        end
        abi = MicrosoftMPIABI

    elseif startswith(MPI_LIBRARY_VERSION_STRING, "Intel")
        impl = IntelMPI

        # TODO: figure out how to parse
        # "Intel(R) MPI Library 2019 Update 4 for Linux* OS"
        if (m = match(r"^Intel\(R\) MPI Library (\d+)", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
            if version > v"2014"
                abi = MPICHABI
            end
        end

    elseif startswith(MPI_LIBRARY_VERSION_STRING, "MVAPICH2")
        impl = MVAPICH
        # "MVAPICH2 Version      :\t%s\n")
        if (m = match(r"^MVAPICH2? Version\s*:\t(\S*)\n", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
            if version > v"2"
                abi = MPICHABI
            end
        end
    end

    if (abienv = get(ENV, "JULIA_MPI_ABI", nothing)) !== nothing
        for inst in instances(MPIABI)
            if String(Symbol(inst)) == abienv
                abi = inst
            end
        end
    end

    return impl, version, abi
end

const MPI_LIBRARY, MPI_LIBRARY_VERSION, MPI_LIBRARY_ABI = identify_implementation()

"""
    MPI_LIBRARY :: MPIImpl

The current MPI implementation: this is determined by

# See also
- [`MPIImpl`](@ref)
"""
MPI_LIBRARY

"""
    MPI_LIBRARY_VERSION :: VersionNumber

The version of the MPI library
"""
MPI_LIBRARY_VERSION

"""
    MPI_LIBRARY_ABI :: MPIABI

The ABI used by the current MPI implementation.

# See also
- [`MPIABI`](@ref)
"""
MPI_LIBRARY_ABI

function Get_version()
    major = Ref{Cint}()
    minor = Ref{Cint}()
    @mpicall ccall((:MPI_Get_version, libmpi), Cint,
                   (Ptr{Cint}, Ptr{Cint}), major, minor)
    VersionNumber(major[], minor[])
end

"""
    MPI_VERSION :: VersionNumber

The supported version of the MPI standard.

# External links
$(_doc_external("MPI_Get_version"))
"""
const MPI_VERSION = Get_version()


text:
module IterativeRefinement
# this file is part of IterativeRefinement.jl, released under the MIT Expat license.

using LinearAlgebra

export rfldiv, equilibrators, condInfest, rfeigen

include("infra.jl")

# Algorithm 3 from
# J.Demmel et al., "Error bounds from extra precise iterative refinement",
# LAPACK Working Note Nr. 165 (2005), also published as
# ACM TOMS, 32, 325 (2006) (henceforth "the paper").

"""
    rfldiv(A,b,f=lu; kwargs...) -> x,bnorm,bcomp,flags

Compute an accurate solution to a linear system ``A x = b`` using
extra-precise iterative refinement, with error bounds.

Returns solution `x`, a normwise relative forward error estimate `bnorm`,
and maximum componentwise relative error estimate `bcomp`.
Specifically,  `bnorm` is an estimate of  ``‖xtrue - x‖ / ‖x‖`` (max norms).
If the problem is so ill-conditioned that a good solution is unrealizable,
`bnorm` and `bcomp` are set to unity (unless `expert`).
`flags` contains convergence diagnostics potentially interesting to
specialists.

# Arguments
- `A`: a matrix,
- `b`: a vector with the same `eltype`,
- `f`: a factorization function such as `lu`.

## Keywords

- `DT`: higher-precision type for refinement; defaults to `widen(eltype(A))`
- `verbosity`: 0(default): quiet, 1: report on iterations, 2: details.
- `equilibrate::Bool`: whether the function should equilibrate `A`
  (default `true`).
- `maxiter`: default 20.
- `tol`: relative tolerance for convergence, in units of `eps(T)`.
- `expert::Bool`: whether to return questionable bounds in extreme cases.
- `κ`: the (max-norm) condition of `A` (see below).
- `F`: a factorization of `A` (see below).

If `A` has already been equilibrated, and a `Factorization` object `F`
and condition estimate `κ` have already been computed, they may be
provided as keyword arguments; no check for consistency is done here.

Uses the algorithm of Demmel et al. ACM TOMS, 32, 325 (2006).
"""
function rfldiv(A::AbstractMatrix{T},
                b::AbstractVecOrMat{T},
                factor = lu;
                DT = widen(T),
                maxiter=20, tol=max(10.0,sqrt(size(A,1))),
                equilibrate = true,
                verbosity = 0,
                ρthresh = 0.5, # "what Wilkinson used"
                expert = false,
                κ = -one(real(T)),
                F::Union{Nothing, Factorization} = nothing
                ) where {T}
    RT = real(T)
    rfldiv_(A,b,lu,DT,RT,maxiter,tol,equilibrate,verbosity,ρthresh,expert,κ,F)
end
function rfldiv_(A::AbstractMatrix{T},
                 b::AbstractVecOrMat{T},
                 factor,
                 ::Type{DT},
                 ::Type{RT},
                 maxiter, tol, equilibrate, verbosity, ρthresh, expert, κ, F
                ) where {T, DT, RT}
    # maxiter is ithresh in paper
    # tol is γ in paper

    m,n = size(A,1), size(A,2)
    if size(b,1) != m
        throw(DimensionMismatch("first dimension of A, $n, does not match that of b, $(size(b,1))"))
    end
    nrhs = size(b,2)

    cvtok = true

    # the use of this variable in closures subverts inference, as of v1.1
    ϵw::RT = RT(2)^(-precision(RT)) # typically eps(T) / 2
    tol1 = 1 / (tol * ϵw) # $1/γϵ_w$ in the paper

    if equilibrate
        Rv,Cv = equilibrators(A)
        cnorm = maximum(abs.(Cv))
        equil = cnorm > 10
    else
        equil = false
    end
    if equil
        (verbosity > 1) && println("equilibrating, ‖C‖=$cnorm")
        C = Diagonal(Cv)
        R = Diagonal(Rv)
        As = R * A * C
    else
        C = I
        As = A
    end

    local Asd
    try
        Asd = DT.(As)
    catch
        cvtok = false
    end
    cvtok || throw(ArgumentError("unable to convert to "
                                 * "designated wide type $DT"))

    if F === nothing
        Fs = factor(As)
    else
        Fs = F
    end
    if κ < 0
        anorm = opnorm(As, Inf)
        κs = condInfest(As,Fs,anorm)
        if verbosity > 1
            equil && print("equilibrated ")
            println("norm: $anorm condition: $κs; compare to $tol1")
        end
    else
        κs = κ
    end

    dzthresh = 1/4 # "works well for binary arithmetic"

    nsingle = 1
    ndouble = 0
    relnormx = relnormz = RT(Inf)
    dxnorm = dznorm = RT(Inf)
    ρmax_x = ρmax_z = zero(RT)
    xstate = :working
    zstate = :unstable
    yscheme = :single
    incrprec = false

    if nrhs > 1
        X = zeros(T,n,nrhs)
        normwisebounds = zeros(RT,nrhs)
        termwisebounds = zeros(RT,nrhs)
        flagss = zeros(Int,3,nrhs)
    end

    function newxstate(state, xnorm, dxnorm, dxprev)
        curnorm = relnormx
        dxratio = dxnorm / dxprev
        dxrel = dxnorm / xnorm
        if (state == :noprogress) && (dxratio <= ρthresh)
            state = :working
        end
        if state == :working
            if dxrel <= ϵw
                # tiny dx, criterion (18) in paper
                state = :converged
                (verbosity > 1) && println("convergence (in norm)")
            elseif dxratio > ρthresh
                if yscheme == :single
                    (verbosity > 1) && println("increasing precision(x)")
                    incrprec = true
                elseif ndouble > 1
                    # lack of progress, criterion (17) in paper
                    state = :noprogress
                    (verbosity > 1) && println("x stalled")
                end
            else
                ρmax_x = max( ρmax_x, dxratio)
            end
            (state != :working) && (curnorm = dxrel)
        end
        state, ρmax_x, curnorm
    end
    function newzstate(state, dznorm, dzprev)
        curnorm = relnormz
        dzratio = dznorm / dzprev
        if (state == :unstable) && (dznorm <= dzthresh)
            state = :working
        end
        if (state == :noprogress) && (dzratio <= ρthresh)
            state = :working
        end
        if state == :working
            if dznorm <= ϵw
                # tiny dz
                state = :converged
                (verbosity > 1) && println("convergence (component-wise)")
            elseif dznorm > dzthresh
                state = :unstable
                relnormz = RT(Inf)
                ρmax_z = zero(RT)
            elseif dzratio > ρthresh
                if yscheme == :single
                    (verbosity > 1) && println("increasing precision(z)")
                    incrprec = true
                elseif ndouble > 1
                    state = :noprogress
                    (verbosity > 1) && println("z stalled")
                end
            else
                ρmax_z = max(ρmax_z, dzratio)
            end
            (state != :working) && (curnorm = dznorm)
        end
        state, ρmax_z, curnorm
    end

    # simple "for" loop w/o scoping
    irhs = 1

    @label NEXTRHS

    if equil
        bs = R * b[:,irhs]
    else
        bs = b[:,irhs]
    end
    bd = DT.(bs)
    y = Fs \ bs

    local yd, xnorm
    for iter=1:maxiter
        # compute residual in appropriate precision
        if yscheme == :single
            r = As * y - bs
            nsingle += 1
        else
            r = T.(Asd * yd - bd)
            ndouble += 1
        end
        # compute correction to y
        dy = Fs \ r
        # check error-related stopping criteria
        xnorm = norm(C*y,Inf)
        dxprev = dxnorm
        dxnorm = norm(C*dy,Inf)
        dzprev = dznorm
        dznorm = maximum( abs.(dy) ./ abs.(y))
        (verbosity > 0) && println("iter $iter |dx|=$dxnorm |dz|=$dznorm")
        ay0,ay1 = extrema(abs.(y))
        if (yscheme == :single) && (κs * ay1 / ay0 >= tol1)
            (verbosity > 1) && println("increasing precision")
            incrprec = true
        end
        xstate, ρmax_x, relnormx = newxstate(xstate, xnorm, dxnorm, dxprev)
        zstate, ρmax_z, relnormz = newzstate(zstate, dznorm, dzprev)
        # the unstable z case is not in the paper but seems
        # necessary to prevent early stalling
        if ((xstate != :working) && !(zstate ∈ [:working,:unstable]))
            break
        end
        if incrprec
# with modified logic above:
#            if yscheme == :double
#                @warn "secondary widening is indicated but not implemented"
#            end
            yscheme = :double
            incrprec = false
            yd = DT.(y)
        end
        # update solution
        if yscheme == :single
            y .-= dy
        else
            yd .-= DT.(dy)
            y = T.(yd)
        end
    end
    if xstate == :working
        relnormx = dxnorm / xnorm
    end
    if zstate == :working
        relnormz = dznorm
    end
    x::Vector{T} = C * y
    min1 = max(10,sqrt(n)) * ϵw # value from paper
    min2 = ϵw
    normwisebound = RT(max( relnormx/(1-ρmax_x), min2))
    termwisebound = RT(max( relnormz/(1-ρmax_z), min1))
    if !expert
        flag = false
        if normwisebound > sqrt(ϵw)
            flag = true
            normwisebound = one(RT)
        end
        if termwisebound > sqrt(ϵw)
            flag = true
            termwisebound = one(RT)
        end
        if flag && (verbosity >= 0)
            @warn "no convergence: result is not meaningful"
        end
    end
    fval = Dict(:converged => 0, :working => 1,
                :noprogress => 2, :unstable => 3)
    flags = [10*fval[xstate]+fval[zstate],nsingle,ndouble]

    # let's see if we can make this type-stable.
    if !(b isa AbstractVector)

        X[:,irhs] .= x
        normwisebounds[irhs] = normwisebound
        termwisebounds[irhs] = termwisebound
        flagss[:,irhs] .= flags

        if irhs == nrhs
            return (X, normwisebounds, termwisebounds, flagss)
        end

    else
        # return (X, normwisebounds, termwisebounds, flagss)
        # take that, you poor confused compiler!
        nb::RT = convert(RT,normwisebound)
        tb::RT = convert(RT,termwisebound)
        ff::Vector{Int} = convert.(Int,flags)
        return x, nb, tb, ff
    end

    # re-initialize state
    nsingle = 1
    ndouble = 0
    relnormx = relnormz = RT(Inf)
    dxnorm = dznorm = RT(Inf)
    ρmax_x = ρmax_z = zero(RT)
    xstate = :working
    zstate = :unstable
    yscheme = :single
    incrprec = false

    irhs += 1
    @goto NEXTRHS
    # no path to this location
end

include("eigen.jl")
end # module


text:
using LabelledArrays

export chuacircuit_nscroll_sine

@inline @inbounds function eom_chuacircuit_nscroll_sine(u, p, t)
    α, β, γ, a, b, c, σx, σy, σz = (p...,)
    x, y, z = (u...,)

    n::Int = c + 1
    if x >= 2*a*c
        fx = (b*pi/2*a)*(x - 2*a*c)
    elseif -2*a*c < x < 2*a*c
        d = ifelse(isodd(n), pi, 0)
        fx = -b*sin((pi*x/2*a) + d)
    elseif x <= -2*a*c
        fx = (b*pi/2*a)*(x + 2*a*c)
    end

    ηx = σx == 0 ? 0 : rand(Normal(0, σx))
    ηy = σy == 0 ? 0 : rand(Normal(0, σy))
    ηz = σz == 0 ? 0 : rand(Normal(0, σz))

    dx = α*(y - fx) + ηx
    dy = x - y + z + ηy
    dz = -β*y - γ*z + ηz
    return SVector{3}(dx, dy, dz)
end

function chuacircuit_nscroll_sine(u₀, α, β, γ, a, b, c::Int, σx, σy, σz)
    p = @LArray [α, β, γ, a, b, c, σx, σy, σz] (:α, :β, :γ, :a, :b, :c, :σx, :σy, :σz)
    ContinuousDynamicalSystem(eom_chuacircuit_nscroll_sine, u₀, p)
end

"""
    chuacircuit_nscroll_sine(;u₀ = [0.0, 0.0, 0.28695],
        α = 10.814, β = 14, γ = 0, a = 1.3, b = 0.11, c = 2,
        σx = 0.0, σy = 0.0, σz = 0.0)

Initialise an adjusted Chua system giving rise to n-scroll attractors [1].

## Equations of motion

The dynamics is generated by the following vector field

```math
\\begin{aligned}
\\dot{x} = \\alpha (y - fx) + \\eta x \\\\
\\dot{y} = x - y + z + \\eta y \\\\
\\dot{z} = -\\beta y - \\gamma z + \\eta z
\\end{aligned}
```

where ``\\eta x``, ``\\eta z``, and ``\\eta z`` are drawn independently from
normal distributions with zero mean and standard deviations `σx`, `σy`
and `σz` at each iteration.

``fx`` is given by the following conditions:

```julia
n::Int = c + 1

if x >= 2*a*c
    fx = (b*pi/2*a)*(x - 2*a*c)
elseif -2*a*c < x < 2*a*c
    d = ifelse(isodd(n), pi, 0)
    fx = -b*sin((pi*x/2*a) + d)
elseif x <= -2*a*c
    fx = (b*pi/2*a)*(x + 2*a*c)
end
```

## References

1. Tang, Wallace KS, et al. "Generation of n-scroll attractors via
    sine function." IEEE Transactions on Circuits and Systems I:
    Fundamental Theory and Applications 48.11 (2001): 1369-1372.
"""
chuacircuit_nscroll_sine(;u₀ = [0.0, 0.0, 0.28695],
        α = 10.814, β = 14, γ = 0, a = 1.3, b = 0.11, c = 2,
        σx = 0.0, σy = 0.0, σz = 0.0) =
    chuacircuit_nscroll_sine(u₀, α, β, γ, a, b, c, σx, σy, σz)


text:
using Documenter, EMST

makedocs(modules=[EMST],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/lithom/EMST.git",
    julia  = "0.6.2",
    osname = "Windows")


text:
export ManualActuator

"""
`AbstractManualActuator`

An abstract interface to actuators with manual input.
"""
abstract type AbstractManualActuator <: AbstractActuator end

mutable struct ManualActuator{MSG} <: AbstractManualActuator
    devname::String
    "Name variable commanded by `ManualActuator`"
    var::String
    "Present value of the actuator"
    val::Float64
    "Minimum value of `val`"
    minval::Float64
    "Maximum value of `val`"
    maxval::Float64
    "Interface to be used to actually wait for user input"
    interface::MSG
    "Time in seconds to wait after user presses ENTER"
    nsec::Float64
end

"""
`TermMSG`

Simple interface using `stdin` to get input from user
"""
struct TermMSG
    msg::String
    confirm::Bool
    TermMSG(msg="", confirm=true) = new(msg, confirm)
end

function (tmsg::TermMSG)(dev::T, x) where {T<:AbstractManualActuator}
    println(tmsg.msg)
    print("Set $(dev.devname) = $x and press ENTER to continue...")
    readline()

    if tmsg.confirm
        print("If ready, press ENTER to start...")
        readline()
    end

end

"""
`ManualActuator(devname, var, val, interf; minval=Inf, maxval=Inf, nsec=0.0)`

`ManualActuator(devname, var, val; minval=Inf, maxval=Inf, nsec=0.0)`

Creates a manual actuator. A manual actuator is an actuator that is not
automated. The operator should manually operate it.

This interface basically explicitly asks the operator to operate the actuator
and waits until the operator says it can go on.

The issue here is the interface. The [`TermMSG`](@ref) uses the terminal to
show information and get input from the operator. Other interfaces such as
different types of GUI could be implemented.

The only requirement is that the object of type `MSG` should be callable with
two arguments:

 * The `ManualActuator` itself
 * The position that the actuator should go to.

This function will get input from the user to do this.

"""
function ManualActuator(devname, var, val, interf::MSG; minval=-Inf,
                        maxval=Inf, nsec=0.0) where {MSG}
    ManualActuator{MSG}(devname, var, val, minval, maxval, interf, nsec)
end

ManualActuator(devname, var, val=0; minval=-Inf, maxval=Inf, nsec=0.0) =
    ManualActuator(devname, var, val, TermMSG(); minval=minval, maxval=maxval, nsec=nsec)

numaxes(dev::ManualActuator) = 1
axesnames(dev::ManualActuator) = [dev.var]

function move(dev::ManualActuator{T}, x) where{T}

    if x < dev.minval || x > dev.maxval
        throw(DomainError(x, "Outside valid range ($(dev.minval), $(dev.maxval))"))
    end
    dev.interface(dev, x)
    dev.val = x
    sleep(dev.nsec)
    return

end

moveto(dev::ManualActuator{T}, x) where{T} = move(dev, x[1])
devposition(dev::ManualActuator) = dev.val


text:
"""
    stairs(x, y; kwargs...)

Description
============

Draws a staircase plot on a new canvas.

The first (optional) vector `x` should contain the horizontal
positions for all the points. The second vector `y` should then
contain the corresponding vertical positions respectively. This
means that the two vectors must be of the same length and
ordering.

Usage
======

    stairs(x, y; style = :post, name = "", title = "", xlabel = "", ylabel = "", labels = true, border = :solid, margin = 3, padding = 1, color = :auto, width = 40, height = 15, xlim = (0, 0), ylim = (0, 0), canvas = BrailleCanvas, grid = true)

Arguments
==========

- **`x`** : The horizontal position for each point.

- **`y`** : The vertical position for each point.

- **`style`** : Specifies where the transition of the stair takes
  plays. Can be either `:pre` or `:post`.

- **`name`** : Annotation of the current drawing to be displayed
  on the right.

$DOC_PLOT_PARAMS

- **`height`** : Number of character rows that should be used
  for plotting.

- **`xlim`** : Plotting range for the x axis.
  `(0, 0)` stands for automatic.

- **`ylim`** : Plotting range for the y axis.
  `(0, 0)` stands for automatic.

- **`canvas`** : The type of canvas that should be used for drawing.

- **`grid`** : If `true`, draws grid-lines at the origin.

Returns
========

A plot object of type `Plot{T<:Canvas}`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)
- Dominique (Github: https://github.com/dpo)

Examples
=========

```julia-repl
julia> stairs([1, 2, 4, 7, 8], [1, 3, 4, 2, 7], style = :post, title = "My Staircase Plot")
                 My Staircase Plot
     ┌────────────────────────────────────────┐
   7 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⡄⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⢸⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠧⠤⠤⠤⠤⠼│
     │⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   1 │⣀⣀⣀⣀⣀⣸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
     └────────────────────────────────────────┘
     1                                        8
```

See also
=========

[`Plot`](@ref), [`scatterplot`](@ref), [`lineplot`](@ref),
[`BrailleCanvas`](@ref), [`BlockCanvas`](@ref),
[`AsciiCanvas`](@ref), [`DotCanvas`](@ref)
"""
function stairs(X::AbstractVector, Y::AbstractVector; style::Symbol = :post, kw...)
    x_vex, y_vex = compute_stair_lines(X, Y, style)
    lineplot(x_vex, y_vex; kw...)
end

function stairs!(
    plot::Plot{<:Canvas},
    X::AbstractVector,
    Y::AbstractVector;
    style::Symbol = :post,
    kw...,
)
    x_vex, y_vex = compute_stair_lines(X, Y, style)
    lineplot!(plot, x_vex, y_vex; kw...)
end

function compute_stair_lines(X::AbstractVector, Y::AbstractVector, style::Symbol)
    if style == :post
        x_vex = similar(X, 2length(X) - 1)
        y_vex = similar(Y, 2length(X) - 1)
        x_vex[1] = X[1]
        y_vex[1] = Y[1]
        o = 0
        for i in 2:(length(X))
            x_vex[i + o] = X[i]
            x_vex[i + o + 1] = X[i]
            y_vex[i + o] = Y[i - 1]
            y_vex[i + o + 1] = Y[i]
            o += 1
        end
        return x_vex, y_vex
    elseif style == :pre
        x_vex = zeros(2length(X) - 1)
        y_vex = zeros(2length(X) - 1)
        x_vex[1] = X[1]
        y_vex[1] = Y[1]
        o = 0
        for i in 2:(length(X))
            x_vex[i + o] = X[i - 1]
            x_vex[i + o + 1] = X[i]
            y_vex[i + o] = Y[i]
            y_vex[i + o + 1] = Y[i]
            o += 1
        end
        return x_vex, y_vex
    end
end


text:
__precompile__()

module SymEngine

import Base: show, convert, real, imag
import Compat: String, unsafe_string, @compat, denominator, numerator, invokelatest

export Basic, symbols, @vars, @funs, SymFunction
export free_symbols, get_args
export ascii_art
export subs, lambdify, N
export series

const deps_file = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
const deps_in_file = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl.in")

noinit = false

if !isfile(deps_file)
    cp(deps_in_file, deps_file)
end

include(deps_file)
include("utils.jl")
const have_mpfr = have_component("mpfr")
const have_mpc = have_component("mpc")
const libversion = get_libversion()

include("types.jl")
include("ctypes.jl")
include("display.jl")
include("mathops.jl")
include("mathfuns.jl")
include("subs.jl")
include("numerics.jl")
include("simplify.jl")
include("calculus.jl")
include("recipes.jl")
include("dense-matrix.jl")

if !noinit
    __init__() = init_constants()
end

end


text:
# initial tests for readTopology
# reads only a tree and with taxon number, not taxon names
# Claudia October 2014
#########################################################

include("../src/types.jl")
include("../src/functions.jl")

using Base.Collections # for updateInCycle with priority queue

# good trees ---------------------------
f = open("prueba_tree.txt","w")
tree = "((1,2),(3,4));"
tree = "((11,22),(33,44));"
tree = "((Ant,Bear),(Cat,Dog));"
tree = "((Ant1,Bear2),(Cat3,Dog4));"
tree = "((1Ant,2Bear),(3Cat,4Dog));"
tree = "((1,2),3,4);"
tree = "(1,2,(3,4));"
tree = "(Ant,Bear,(Cat,Dog));"
tree = "(A,B,(C,D));"
tree = "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"
tree = "((((((1,2),3),4),5),6),7,8);" # yeast data tree0
write(f,tree)
close(f)

net = readTopology("prueba_tree.txt");
printEdges(net)
printNodes(net)
net.names

cleanAfterRead!(net)

# bad trees --------------------------
f = open("prueba_tree.txt","w")
tree = "(((1,2),(3,4)));" # extra parenthesis
tree = "((1,2,(3,4)));" # extra parenthesis
tree = "(((1,2),3,4));" # extra parenthesis
tree = "((1),2,3,4,5);" # not tree
tree = "((1,2),3,4,5);" # polytomy
tree = "((1,*),(3,4));" # not letter/number taxon name
tree = "(1,2::,(3,4));" #double :
tree = "(1,2:::,(3,4));" #triple :
tree = "(1,2:0.2:,(3,4));" #no 2nd :
tree = "(1,2:0.2:0.2:,(3,4));" #no 3rd :
tree = "(1,2:,(3,4));" #no 1st:
tree = "((1,2),(3,4))" # no ;
write(f,tree)
close(f)

net = readTopology("prueba_tree.txt")
printEdges(net)
printNodes(net)
net.names

cleanAfterRead!(net)

# trees with additional info -----------------
f = open("prueba_tree.txt","w")
tree = "((1:1.2,2:0.3),3:1.8,4);"
tree = "((1:1.2,2:0.3):1.2:0.6:0.3,3:1.8,4);"
tree = "((1:1.2,2:0.3):1.2::0.6,3:1.8,4);" # no bootstrap
tree = "((1:1.2,2:0.3):::0.6,3:1.8,4);" # no length nor bootstrap
tree = "((1:1.2,2:0.3):1.2:2.5:1.6,3:1.8,4);" # wrong value gamma
tree = "((1:1.2,2:0.3):1.2:2.5:,3:1.8,4);" # missing value gamma
write(f,tree)
close(f)

net = readTopology("prueba_tree.txt")
printEdges(net)
printNodes(net)
net.names

# beginning of networks ---------------------------
include("../src/types.jl")
include("../src/functions.jl")

using Base.Collections # for updateInCycle with priority queue

f = open("prueba_tree.txt","w")
tree = "(((3,4)Z#H1,1),(Z#H1,2));" # expand child, no gammas
tree = "((Z#H1,2),((3,4)Z#H1,1));" # expand child, no gammas, leaf read first
tree = "(((3,4)Z#H1:::0.9,1),(Z#H1:::0.2,2));" # gammas do not sum up to 1
tree = "(((3,4)Z#H1:5.0::0.9,1),(Z#H1:::0.2,2));"
tree = "(1,2,(3,4)A);"
tree = "(1,2,(3,4)A:0.8);"
tree = "((1,2),#H1,4);" #only one hybrid
tree = "((1,2),#H1,(3,#H2));" #two hybrids but different
tree = "((1,2)#H1,(3,4)#H1);" #error: both hybrid have childre
tree = "(((1,#H1),2),#H1,3);" # both hybrid leaves
tree = "(((1,4)#H1,2),#H1,(5,6)#H1);" # hybrid polytomy
write(f,tree)
close(f)

net = readTopology("prueba_tree.txt");
printEdges(net)
printNodes(net)
net.names
net.numHybrids

updateAllReadTopology!(net)

readTopologyUpdate("prueba_tree.txt")

#======================================================

include("../src/types.jl")
include("../src/functions.jl")

using Base.Collections # for updateInCycle with priority queue

# good trees ---------------------------
f = open("prueba_tree.txt","w")
tree = "((1,2),(3,4));"
tree = "((11,22),(33,44));"
tree = "((Ant,Bear),(Cat,Dog));"
tree = "((Ant1,Bear2),(Cat3,Dog4));"
tree = "((1Ant,2Bear),(3Cat,4Dog));"
tree = "((1,2),3,4);"
tree = "(1,2,(3,4));"
tree = "(Ant,Bear,(Cat,Dog));"
tree = "(A,B,(C,D));"
tree = "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"
tree = "((((((1,2),3),4),5),6),7,8);" # yeast data tree0
tree = "(Scer,((Smik,(Skud,Sbay)),Spar));" # yeast data astral output
tree = "(1,2,3,4,5);"
write(f,tree)
close(f)

net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)
printNodes(net)
net.names

net2 = readTopologyUpdate(tree);
printEdges(net2)
printNodes(net2)
net2.names


text:
# Original implementations by @nantonel
# https://github.com/maxmouchet/HMMBase.jl/pull/6

# *log! methods use the samples log-likelihood instead of the likelihood.

# In-place forward pass, where α and c are allocated beforehand.
function forward!(α::AbstractMatrix, c::AbstractVector, a::AbstractVector, A::AbstractMatrix, L::AbstractMatrix)
    @argcheck size(α, 1) == size(L, 1) == size(c, 1)
    @argcheck size(α, 2) == size(L, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(L)
    (T == 0) && return

    fill!(α, 0.0)
    fill!(c, 0.0)

    for j in OneTo(K)
        α[1,j] = a[j] * L[1,j]
        c[1] += α[1,j]
    end

    for j in OneTo(K)
        α[1,j] /= c[1]
    end

    @inbounds for t in 2:T
        for j in OneTo(K)
            for i in OneTo(K)
                α[t,j] += α[t-1,i] * A[i,j]
            end
            α[t,j] *= L[t,j]
            c[t] += α[t,j]
        end

        for j in OneTo(K)
            α[t,j] /= c[t]
        end
    end
end

# In-place backward pass, where β and c are allocated beforehand.
function backward!(β::AbstractMatrix, c::AbstractVector, a::AbstractVector, A::AbstractMatrix, L::AbstractMatrix)
    @argcheck size(β, 1) == size(L, 1) == size(c, 1)
    @argcheck size(β, 2) == size(L, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(L)
    (T == 0) && return

    fill!(β, 0.0)
    fill!(c, 0.0)

    for j in OneTo(K)
        β[end,j] = 1.0
    end

    @inbounds for t in T-1:-1:1
        for j in OneTo(K)
            for i in OneTo(K)
                β[t,j] += β[t+1,i] * A[j,i] * L[t+1,i]
            end
            c[t+1] += β[t,j]
        end

        for j in OneTo(K)
            β[t,j] /= c[t+1]
        end
    end

    for j in OneTo(K)
        c[1] += a[j] * L[1,j] * β[1,j]
    end
end

# In-place forward pass, where α and c are allocated beforehand.
function forwardlog!(α::AbstractMatrix, c::AbstractVector, a::AbstractVector, A::AbstractMatrix, LL::AbstractMatrix)
    @argcheck size(α, 1) == size(LL, 1) == size(c, 1)
    @argcheck size(α, 2) == size(LL, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(LL)
    (T == 0) && return

    fill!(α, 0.0)
    fill!(c, 0.0)

    m = vec_maximum(view(LL, 1, :))

    for j in OneTo(K)
        α[1,j] = a[j] * exp(LL[1,j] - m)
        c[1] += α[1,j]
    end

    for j in OneTo(K)
        α[1,j] /= c[1]
    end

    c[1] *= exp(m) + eps()

    @inbounds for t in 2:T
        m = vec_maximum(view(LL, t, :))

        for j in OneTo(K)
            for i in OneTo(K)
                α[t,j] += α[t-1,i] * A[i,j]
            end
            α[t,j] *= exp(LL[t,j] - m)
            c[t] += α[t,j]
        end

        for j in OneTo(K)
            α[t,j] /= c[t]
        end

        c[t] *= exp(m) + eps()
    end
end

# In-place backward pass, where β and c are allocated beforehand.
function backwardlog!(β::AbstractMatrix, c::AbstractVector, a::AbstractVector, A::AbstractMatrix, LL::AbstractMatrix)
    @argcheck size(β, 1) == size(LL, 1) == size(c, 1)
    @argcheck size(β, 2) == size(LL, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(LL)
    L = zeros(K)
    (T == 0) && return

    fill!(β, 0.0)
    fill!(c, 0.0)

    for j in OneTo(K)
        β[end,j] = 1.0
    end

    @inbounds for t in T-1:-1:1
        m = vec_maximum(view(LL, t+1, :))

        for i in OneTo(K)
            L[i] = exp(LL[t+1,i] - m)
        end

        for j in OneTo(K)
            for i in OneTo(K)
                β[t,j] += β[t+1,i] * A[j,i] * L[i]
            end
            c[t+1] += β[t,j]
        end

        for j in OneTo(K)
            β[t,j] /= c[t+1]
        end

        c[t+1] *= exp(m) + eps()
    end

    m = vec_maximum(view(LL, 1,:))

    for j in OneTo(K)
        c[1] += a[j] * exp(LL[1,j] - m) * β[1,j]
    end

    c[1] *= exp(m) + eps();
end

# In-place posterior computation, where γ is allocated beforehand.
function posteriors!(γ::AbstractMatrix, α::AbstractMatrix, β::AbstractMatrix)
    @argcheck size(γ) == size(α) == size(β)
    T, K = size(α)
    for t in OneTo(T)
        c = 0.0
        for i = OneTo(K)
            γ[t,i] = α[t,i] * β[t,i]
            c += γ[t,i]
        end

        for i in OneTo(K)
            γ[t,i] /= c
        end
    end
end

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
function convexHull( points::Array{Tuple{T,T,T},1} ) where {T<:Real}

    @assert length(points) >= 4
    # length(points) for now, to have a static size

    # Each face consistes of 3 Vertices and 3 edges.
	# The edges of a face, face[idx], are given by face[1]-face[2], face[2]-face[3], face[3]-face[1];
    faces = Array{Tuple{UInt16,UInt16,UInt16},1}(undef,length(points)*2);

    # Each face is in contact with 3 more faces
    neigh = Array{Tuple{UInt16,UInt16,UInt16},1}(undef,length(points)*2);

    # One normal for each face
    norms = Array{Tuple{Float32,Float32,Float32},1}(undef,length(points)*2);

    # Intializing first 4 faces and neighbouring info NOTE: need to assert that 4 first points are not colinear/coplanar
    faces[1] = ( 1, 2, 3 ); neigh[1] = ( 4, 2, 3 ); # edge 1-2 is shared with faces[4], edge 2-3 with faces[2], and edge 3-1 with faces[3]
    faces[2] = ( 2, 3, 4 ); neigh[2] = ( 1, 3, 4 ); # edge 2-3 is shared with faces[1], edge 3-4 with faces[3], and edge 4-2 with faces[4]
    faces[3] = ( 3, 4, 1 ); neigh[3] = ( 2, 4, 1 ); # edge 3-4 is shared with faces[2], edge 4-1 with faces[4], and edge 1-3 with faces[1]
    faces[4] = ( 4, 1, 2 ); neigh[4] = ( 3, 1, 2 ); # edge 4-1 is shared with faces[3], edge 1-2 with faces[1], and edge 2-4 with faces[2]

    # Faces will be modified and added. When adding new faces, "last" keeps track of the index of the last face.
    lastFace = 4;

	center = ( points[1] .+ points[2] .+ points[3] .+ points[4] ) ./ 4;

    areCoplanar( points[1], points[2], points[3], points[4] )  && (println("initial are coplanar"))

    for fidx in 1:4
		norms[fidx] = faceNormalCenter( points, faces[fidx], center )
    end

	# Iteratively modifying the current convex Hull

	for pidx in 5:length(points);

		@inbounds new_point = points[pidx];

		# 1-. Checking faces visible to new_point.
		visibleFaces = Array{UInt16,1}(undef,lastFace);
		numVisibleFaces = UInt16(0);

		for fidx in 1:lastFace;

			areCoplanar( points[faces[fidx][1]], points[faces[fidx][2]], points[faces[fidx][3]], new_point ) && (println("middle are coplanar"))

			@inbounds vp = new_point .- points[faces[fidx][3]];
			@inbounds tf =  sign( sum( vp .* norms[ fidx ] ) );
			if tf == 1
				numVisibleFaces += UInt16(1);
				visibleFaces[numVisibleFaces] = fidx;
			end
		end

		lastFace = addNewFaces!( points, faces, neigh, norms, center, lastFace,
						         pidx, visibleFaces, numVisibleFaces );

	end # each input point

    return faces[1:lastFace], neigh[1:lastFace], norms[1:lastFace]
end

"""
	needs to be tested --> tested in TEST 1
"""
@inline function faceNormalCenter( points, face, ctr )
    v =  points[face[2]] .- points[face[1]];
    w =  points[face[3]] .- points[face[1]];
	n = ( v[2]*w[3] - v[3]*w[2], v[3]*w[1] - v[1]*w[3], v[1]*w[2] - v[2]*w[1] )
	mag = sqrt( sum( n .* n ) );
	n = n ./ mag
	p = ctr .- points[face[1]];
	return ( sign( sum( p .* n ) ) == -1 ) ? n : -1 .* n;
end

"""
	needs to be tested
"""
@inline function addNewFaces!( points,
                               faces::Array{Tuple{UInt16,UInt16,UInt16},1},
                               neigh::Array{Tuple{UInt16,UInt16,UInt16},1},
							   norms::Array{Tuple{Float32,Float32,Float32},1},
							   center,
							   lastFace,
							   newpointindex,
							   visibleFaces::Array{UInt16,1},
                               numVisibleFaces
                             )
	"""
		1-. Collecting all edges from VisibleFaces in "boundary" and labelling those that are Repeated (not unique)
	"""
	boundary   = Array{Tuple{UInt16,UInt16},1}(undef, numVisibleFaces*3);
	isRepeated = zeros( Bool, numVisibleFaces*3 );

	lidx = 1; # linear index, it goes from 1 to numVisibleFaces*3
	nrep = 0; # number of repeats

	for face in faces[ visibleFaces[1:numVisibleFaces] ]
		@inbounds edges = [(face[1],face[2]), (face[2],face[3]), (face[3],face[1])];
		for (p1, p2) in edges
			@inbounds boundary[lidx] = ( p1, p2 )
			for bidx in 1:lidx-1
				@inbounds boundary_edge = boundary[bidx]
				if (p1,p2) == boundary_edge || (p2,p1) == boundary_edge # && !isRepeated[bidx] # each edge can only be shared by two faces, so this should not be necessary
					isRepeated[lidx] = true;
					isRepeated[bidx] = true;
				end
			end # each previous boundary --> check that the current edge has not already been added
		lidx += 1
		end # each of the 3 edges --> add it to "boundary" and check if it is not repeated
	end # each visible face --> go through its 3 edges

	"""
		2-. From each unique edge in the boundary -> build a new face with the correct neighbouring data
	"""
	# NOTE: lidx is added 1 too much, thus "lidx - 1" below
	nrep = sum( isRepeated ); # nrep += 2 under isRepeated[bidx]
	numNewFaces = lidx - 1 - nrep;

	( numNewFaces < numVisibleFaces ) && ( println( "I Predict shit " ) );

	newfaces = Array{Tuple{UInt16,UInt16,UInt16},1}(undef,numNewFaces);
	newneigh = Array{Tuple{UInt16,UInt16,UInt16},1}(undef,numNewFaces);

	lidx = 1; # linear index, from 1 to numVisibleFaces*3
	nidx = 1; # new face index, from 1 to numNewFaces

	for fidx in visibleFaces[1:numVisibleFaces]
		@inbounds face  = faces[fidx];
		@inbounds edges = [(face[1],face[2]), (face[2],face[3]), (face[3],face[1])];
		for eidx in 1:3

			@inbounds edge = edges[eidx];

			if !isRepeated[lidx]
				# here
				newfaces[nidx] = ( edge[1], edge[2], newpointindex );
				newneigh[nidx] = ( neigh[fidx][eidx], UInt16(0), UInt16(0) );

				newfaceposition = ( nidx <= numVisibleFaces ) ? visibleFaces[nidx] : lastFace + nidx - numVisibleFaces;

				for oidx in 1:nidx-1

					otherface = newfaces[oidx];
					otherfaceposition = ( oidx <= numVisibleFaces ) ? visibleFaces[oidx] : lastFace + oidx - numVisibleFaces;

					# I am sure this can be done smarter
					any( edge .== otherface[1] ) && ( newneigh[ oidx ] = ( newneigh[oidx][1], newneigh[oidx][2],  newfaceposition ) )
					any( edge .== otherface[2] ) && ( newneigh[ oidx ] = ( newneigh[oidx][1],  newfaceposition, newneigh[oidx][3] ) )

				    any( otherface .== edge[1] ) && ( newneigh[nidx] = ( newneigh[nidx][1], newneigh[nidx][2], otherfaceposition ) )
					any( otherface .== edge[2] ) && ( newneigh[nidx] = ( newneigh[nidx][1], otherfaceposition, newneigh[nidx][3] ) )

				end # for other new faces --> add neighbouring data

				nidx += 1;

			end	# for each unique edge --> build a face between that edge and the newpoint

		lidx += 1

		end # for each edge --> check that it is unique
	end # for each visible Face --> go over its 3 edges

	"""
		Either replacing or pushing the new faces with their neighbouring data to "faces" and "neigh"
	"""
	for nidx in 1:numNewFaces
		newfaceposition = ( nidx <= numVisibleFaces ) ? visibleFaces[nidx] : lastFace + nidx - numVisibleFaces;
		# here
		faces[ newfaceposition ] = newfaces[ nidx ]
		neigh[ newfaceposition ] = newneigh[ nidx ]
		norms[ newfaceposition ] = faceNormalCenter( points, newfaces[ nidx ], center )
	end

	if ( numNewFaces < numVisibleFaces )
		tobeDeleted = numVisibleFaces - numNewFaces
		for start in visibleFaces[1+numVisibleFaces-tobeDeleted:numVisibleFaces]
			for fidx in start:lastFace
				faces[fidx] = faces[fidx+1]
			end
		end

		lastFace -= tobeDeleted
		# I still need to fix neighbourhoods
	end

	# C
	#println( lastFace, "  ", numVisibleFaces, "  ", nrep, "  ", numNewFaces, "  ", lastFace + numNewFaces - numVisibleFaces )
	# return new lastFace
	change = ( numNewFaces < numVisibleFaces ) ? 0 : numNewFaces - numVisibleFaces;
	return lastFace + change
end

function areCoplanar( x, y, z, p )

    a1 = x[2] - x[1] ;
    b1 = y[2] - y[1] ;
    c1 = z[2] - z[1] ;
    a2 = x[3] - x[1] ;
    b2 = y[3] - y[1] ;
    c2 = z[3] - z[1] ;
    a = b1 * c2 - b2 * c1 ;
    b = a2 * c1 - a1 * c2 ;
    c = a1 * b2 - b1 * a2 ;
    d = - a * x[1] - b * y[1] - c * z[1];

    return a * p[1] + b * p[2] + c * p[3] + d == 0
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

julia_version = v"1.6.0"

name = "XyceWrapper"
version = v"0.2.0"

# Collection of sources required to complete build
sources = [
    DirectorySource("./src"),
]

# Bash recipe for building across all platforms
script = raw"""
mkdir -p build
cd build
cmake \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DJulia_PREFIX=${prefix} \
    -DCMAKE_BUILD_TYPE=Release \
    ..
cmake --build . --config Release --target install -- -j${nproc}
install_license /usr/share/licenses/MIT
"""

include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libxycelib", :xycelib),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Xyce_jll"; compat="^7.4.0"),
    Dependency("libcxxwrap_julia_jll"),
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("iPredict.jl/src/iPredict.jl")
using .iPredict
using DataFrames
using Plots, Plots.Measures
using JLD2
using DelimitedFiles
using StatsBase

#Initialise variables
GeneIsoformRelationships = "data/gencode_human_gene_isoform_relationships.txt"
NumSimulations = 100
NumIsoformsToSimulate = 4
filteringThreshold = 4
alpha_vals = [0.45, 0.7, 0.95, 1.2, 1.45]
beta_vals = reverse(alpha_vals)
scnorm_counts = ["H1_24_scnorm_counts.txt"]#,
"H1_96_scnorm_counts.txt",
"H9_24_scnorm_counts.txt",
"H9_96_scnorm_counts.txt"]
pFP = 0.01
pFN = 0.04
isoChoiceModels = ["Random", "UniformObserved", "CellVariable"]

for countsMatrix in scnorm_counts

    output = "figures/figure2/" * countsMatrix
    data_out = "figures/figure2_data/" * countsMatrix

    for model in isoChoiceModels

        for i in 1:length(alpha_vals)

            output = "figures/figure2/" * model * string(alpha_vals[i]) * "_" *
            string(beta_vals[i]) * "_" * countsMatrix
            data_out = "figures/figure2_data/" * model * string(alpha_vals[i]) * "_" *
            string(beta_vals[i]) * "_" * countsMatrix

            out = globalPredict("data/" * countsMatrix, GeneIsoformRelationships,
            NumSimulations, output, NumIsoformsToSimulate, filteringThreshold,
            false, alpha_vals[i], beta_vals[i], pFN, pFP, "none", model,
            data_out)
        end
    end
end

function make_hist(path, titleArg)
    Arr = readdlm(path)

    if (titleArg == "3 Isoforms" &&
        path == "figures/figure2_data/H1_96_scnorm_counts.txt_3_isoforms_array.txt")
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1.2),
        ylims = (0, 5),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)

    elseif titleArg == "Real" #|| titleArg == "alpha=0.45, beta=1.45"
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1.2),
        ylims = (0, 5),
        normalize=:pdf,
        xlabel = "Mean No. Isoforms\nper Gene per Cell",
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1.2),
        ylims = (0, 5),
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)
    end
    vline!([mean(Arr)], color = :black, linewidth = 8)
    return p1
end

v0 =  make_hist("figures/figure2_data/Random0.45_1.45_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.45, beta=1.45")
v1 =  make_hist("figures/figure2_data/Random0.7_1.2_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.7, beta=1.2")
v2 =  make_hist("figures/figure2_data/Random0.95_0.95_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.95, beta=0.95")
v3 =  make_hist("figures/figure2_data/Random1.2_0.7_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=1.2, beta=0.7")
v4 =  make_hist("figures/figure2_data/Random1.45_0.45_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=1.45, beta=0.45")

pD = plot(v4, v3, v2, v1, v0,
       layout = (5,1), legend = false, size = (250*8, 600*8),
       axis = (font(180)), titlefontsize=160, margin = 30mm)

savefig(pD, "figures/figure2/Random_figure2D.png")

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

v0 =  make_hist("figures/figure2_data/UniformObserved0.45_1.45_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.45, beta=1.45")
v1 =  make_hist("figures/figure2_data/UniformObserved0.7_1.2_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.7, beta=1.2")
v2 =  make_hist("figures/figure2_data/UniformObserved0.95_0.95_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.95, beta=0.95")
v3 =  make_hist("figures/figure2_data/UniformObserved1.2_0.7_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=1.2, beta=0.7")
v4 =  make_hist("figures/figure2_data/UniformObserved1.45_0.45_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=1.45, beta=0.45")

pD = plot(v4, v3, v2, v1, v0,
       layout = (5,1), legend = false, size = (250*8, 600*8),
       axis = (font(180)), titlefontsize=160, margin = 30mm)

savefig(pD, "figures/figure2/UniformObserved_figure2D.png")

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

v0 =  make_hist("figures/figure2_data/CellVariable0.45_1.45_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.45, beta=1.45")
v1 =  make_hist("figures/figure2_data/CellVariable0.7_1.2_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.7, beta=1.2")
v2 =  make_hist("figures/figure2_data/CellVariable0.95_0.95_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.95, beta=0.95")
v3 =  make_hist("figures/figure2_data/CellVariable1.2_0.7_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=1.2, beta=0.7")
v4 =  make_hist("figures/figure2_data/CellVariable1.45_0.45_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=1.45, beta=0.45")

pD = plot(v4, v3, v2, v1, v0,
       layout = (5,1), legend = false, size = (250*8, 600*8),
       axis = (font(180)), titlefontsize=160, margin = 30mm)

savefig(pD, "figures/figure2/CellVariable_figure2D.png")

####################################################
# OVERLAP
####################################################

function make_hist(path, titleArg)
    Arr = readdlm(path)

    if (titleArg == "3 Isoforms" &&
        path == "figures/figure2_data/H1_96_scnorm_counts.txt_3_isoforms_array.txt")
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        color = palette(:default)[2],
        xlims = (0, 1.2),
        ylims = (0, 5),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)

    elseif titleArg == "Real" #|| titleArg == "alpha=0.45, beta=1.45"
        p1 = histogram(Arr,
        title = titleArg,
        color = palette(:default)[2],
        xlims = (0, 1.2),
        ylims = (0, 5),
        normalize=:pdf,
        xlabel = "Mean No. Isoforms\nper Gene per Cell",
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        color = palette(:default)[2],
        xlims = (0, 1.2),
        ylims = (0, 5),
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)
    end
    vline!([mean(Arr)], color = :black, linewidth = 8)
    return p1
end

v0 =  make_hist("figures/figure2_data/Random0.45_1.45_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.45, beta=1.45")
v1 =  make_hist("figures/figure2_data/Random0.7_1.2_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.7, beta=1.2")
v2 =  make_hist("figures/figure2_data/Random0.95_0.95_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.95, beta=0.95")
v3 =  make_hist("figures/figure2_data/Random1.2_0.7_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=1.2, beta=0.7")
v4 =  make_hist("figures/figure2_data/Random1.45_0.45_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=1.45, beta=0.45")

pD = plot(v4, v3, v2, v1, v0,
       layout = (5,1), legend = false, size = (250*8, 600*8),
       axis = (font(180)), titlefontsize=160, margin = 30mm)

savefig(pD, "figures/figure2/Random_figure2D_overlap.png")

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

v0 =  make_hist("figures/figure2_data/UniformObserved0.45_1.45_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.45, beta=1.45")
v1 =  make_hist("figures/figure2_data/UniformObserved0.7_1.2_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.7, beta=1.2")
v2 =  make_hist("figures/figure2_data/UniformObserved0.95_0.95_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.95, beta=0.95")
v3 =  make_hist("figures/figure2_data/UniformObserved1.2_0.7_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=1.2, beta=0.7")
v4 =  make_hist("figures/figure2_data/UniformObserved1.45_0.45_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=1.45, beta=0.45")

pD = plot(v4, v3, v2, v1, v0,
       layout = (5,1), legend = false, size = (250*8, 600*8),
       axis = (font(180)), titlefontsize=160, margin = 30mm)

savefig(pD, "figures/figure2/UniformObserved_figure2D_overlap.png")

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

v0 =  make_hist("figures/figure2_data/CellVariable0.45_1.45_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.45, beta=1.45")
v1 =  make_hist("figures/figure2_data/CellVariable0.7_1.2_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.7, beta=1.2")
v2 =  make_hist("figures/figure2_data/CellVariable0.95_0.95_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.95, beta=0.95")
v3 =  make_hist("figures/figure2_data/CellVariable1.2_0.7_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=1.2, beta=0.7")
v4 =  make_hist("figures/figure2_data/CellVariable1.45_0.45_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=1.45, beta=0.45")

pD = plot(v4, v3, v2, v1, v0,
       layout = (5,1), legend = false, size = (250*8, 600*8),
       axis = (font(180)), titlefontsize=160, margin = 30mm)

savefig(pD, "figures/figure2/CellVariable_figure2D_overlap.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
using ContinuedFractions, Test

@testset "Continued Fraction" begin
  @test sqrt(2) ≈ continuedfraction(i -> i == 0 ? 1 : 2, x->1, 20)
  # https://oeis.org/A003417
  arrayexp1 = [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1]
  fexp1(i) = arrayexp1[i + 1]
  @test exp(1) ≈ continuedfraction(fexp1, x->1, length(arrayexp1) - 1)
  # https://oeis.org/A001203
  arrayπ = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3]
  fπ(i) = arrayπ[i + 1]
  @test Float64(π) ≈ continuedfraction(fπ, x->1, length(arrayπ) - 1)

  for i ∈ 1:10
    atol, rtol =10.0^-rand(1:15), 10.0^-rand(1:15)
    answer = 1.61803398874989484820
    output = continuedfraction(i->1, x->1; atol=atol, rtol=rtol)[1]
    @test isapprox(answer, output, atol=atol, rtol=rtol)
    answer = sqrt(2)
    output = continuedfraction(i->i==0 ? 1 : 2, x->1; atol=atol, rtol=rtol)[1]
    @test isapprox(answer, output, atol=atol, rtol=rtol)
  end

end