text:
# Types for Julia implementation of pseudolikelihood estimation (Stage2)
# Claudia August 2014
# in julia: include("types.jl")
#
# Types: Edge, Node, HybridNetwork
################################################################

# procedure to create hybrid network:
# 1) create edges defined as hybrid or not
# 2) create nodes with such edges
# 3) setNode! to add nodes into edges
# 4) create hybrid network
# 5) updateGammaz!

# -------------- EDGE -------------------------#
"""
 `ANode`

Abstract node. An object of type [`Edge`](@ref) has a `node` attribute,
which is an vector of 2 `ANode` objects.
The object of type [`Node`](@ref) is an `ANode`, and has an `edge` attribute,
which is vector of `Edge` objects.
"""
abstract type ANode end

"""
    Edge(number)

Data structure for an edge and its various attributes. Most notably:
- `number` (integer): serves as unique identifier;
  remains unchanged when the network is modified,
  with a nearest neighbor interchange for example
- `node`: a vector of [`Node`]s, normally just 2 of them
- `isChild1` (boolean): `true` if `node[1]` is the child node of the edge,
  false if `node[1]` is the parent node of the edge
- `length`: branch length
- `hybrid` (boolean): whether the edge is a tree edge or a hybrid edge
  (in which case `isChild1` is important, even if the network is semi-directed)
- `gamma`: proportion of genetic material inherited by the child node via the edge;
  1.0 for a tree edge
- `isMajor` (boolean): whether the edge is the major path to the child node;
  `true` for tree edges, since a tree edge is the only path to its child node;
  normally true if `gamma>0.5`.

and other fields, used very internally
"""
mutable struct Edge
    number::Int
    length::Float64 #default 1.0
    hybrid::Bool
    y::Float64 # exp(-t), cannot set in constructor for congruence
    z::Float64 # 1-y , cannot set in constructor for congruence
    gamma::Float64 # set to 1.0 for tree edges, hybrid?gamma:1.0
    node::Array{ANode,1}
    isChild1::Bool # used for hybrid edges to set the direction (default true)
    isMajor::Bool  # major edge treated as tree edge for network traversal
                   # true if gamma>.5, or if it is the original tree edge
    inCycle::Int # = Hybrid node number if this edge is part of a cycle created by such hybrid node
    # -1 if not part of cycle. used to add new hybrid edge. updated after edge is part of a network
    containRoot::Bool # true if this edge can contain a root given the direction of hybrid edges
                      # used to add new hybrid edge. updated after edge is part of a network
    istIdentifiable::Bool # true if the parameter t (length) for this edge is identifiable as part of a network
                          # updated after part of a network
    fromBadDiamondI::Bool # true if the edge came from deleting a bad diamond I hybridization
    # inner constructors: ensure congruence among (length, y, z) and (gamma, hybrid, isMajor), and size(node)=2
    Edge(number::Int) = new(number,1.0,false,exp(-1.0),1-exp(-1.0),1.,[],true,true,-1,true,true,false)
    Edge(number::Int, length::Float64) = new(number,length,false,exp(-length),1-exp(-length),1.,[],true,true,-1,true,true,false)
    Edge(number::Int, length::Float64,hybrid::Bool,gamma::Float64) =
        new(number,length,hybrid,exp(-length),1-exp(-length), hybrid ? gamma : 1.,[],true, !hybrid || gamma>0.5 ,-1,!hybrid,true,false)
    Edge(number::Int, length::Float64,hybrid::Bool,gamma::Float64,isMajor::Bool) =
        new(number,length,hybrid,exp(-length),1-exp(-length), hybrid ? gamma : 1.,[],true,isMajor,-1,!hybrid,true,false)
    function Edge(number::Int, length::Float64,hybrid::Bool,gamma::Float64,node::Array{ANode,1})
        size(node,1) != 2 ?
        error("vector of nodes must have exactly 2 values") :
        new(number,length,hybrid,exp(-length),1-exp(-length),
            hybrid ? gamma : 1., node,true, !hybrid || gamma>0.5,
            -1,!hybrid,true,false)
    end
    function Edge(number::Int, length::Float64,hybrid::Bool,gamma::Float64,node::Array{ANode,1},isChild1::Bool, inCycle::Int, containRoot::Bool, istIdentifiable::Bool)
        size(node,1) != 2 ?
        error("vector of nodes must have exactly 2 values") :
        new(number,length,hybrid,exp(-length),1-exp(-length),
            hybrid ? gamma : 1., node,isChild1, !hybrid || gamma>0.5,
            inCycle,containRoot,istIdentifiable,false)
    end
end

# warning: gammaz, inCycle, isBadTriangle/Diamond updated until the node is part of a network
"""
    Node(number, leaf)

Data structure for an edge and its various attributes. Most notably:

- `number` (integer): serves as unique identifier;
  remains unchanged when the network is modified,
  with a nearest neighbor interchange for example
- `leaf` (boolean): whether the node is a leaf (with data typically) or an
  internal node (no data typically)
- `name` (string): taxon name for leaves; internal node may or may not have a name
- `edge`: vector of [`Edge`]s that the node is attached to;
  1 if the node is a leaf, 2 if the node is the root, 3 otherwise, and
  potentially more if the node has a polytomy
- `hybrid` (boolean): whether the node is a hybrid node (with 2 or more parents)
  or a tree node (with a single parent)

Other more internal attributes include:

- `isBadDiamondI` and `isBadDiamondII` (booleans): whether the node is a
  hybrid node where the reticulation forms a cycle of 4 nodes (diamond),
  and where both parents of the hybrid nodes are connected to a leaf.
  In a bad diamond of type I, the hybrid node itself is also connected
  to a leaf but the common neighbor of the 2 hybrid's parents is not connected
  to a leaf.
  In a bad diamond of type II, the hybrid node has an internal node as child,
  and the common neighbor of the 2 hybrid's parents is connected to a leaf.
- `isBadTriangle`, `isVeryBadTriangle` and `isExtBadTriangle` (booleans):
  true if the reticulation forms a cycle of 3 nodes (triangle) and
  depending on the number of leaves attached these 3 nodes. The triangle means
  that the 2 parents of the hybrid node are directly related:
  one is the child of the other. `isBadTriangle` is true if the triangle is
  "good", as per Solís-Lemus & Ané (2016), that is, if all 3 nodes in the cycle
  are not connected to any leaves (the reticulation is detectable from quartet
  concordance factors, even though all branch lengths are not identifiable).
  `isVeryBadTriangle` is true if 2 (or all) of the 3 nodes are connected to a
  leaf, in which case the reticulation is undetectable from unrooted gene tree
  topologies (thus it's best to exclude these reticulations from a search).
  `isBadTriangle` is true if exactly 1 of the 3 nodes is connected to a leaf.

For details see Solís-Lemus & Ané (2016, doi:10.1371/journal.pgen.1005896)
"""
mutable struct Node <: ANode
    number::Int
    leaf::Bool
    hybrid::Bool
    gammaz::Float64  # notes file for explanation. gammaz if tree node, gamma2z if hybrid node.
                     # updated after node is part of network with updateGammaz!
    edge::Array{Edge,1}
    hasHybEdge::Bool #is there a hybrid edge in edge? only needed when hybrid=false (tree node)
    isBadDiamondI::Bool # for hybrid node, is it bad diamond case I, update in updateGammaz!
    isBadDiamondII::Bool # for hybrid node, is it bad diamond case II, update in updateGammaz!
    isExtBadTriangle::Bool # for hybrid node, is it extremely bad triangle, udpate in updateGammaz!
    isVeryBadTriangle::Bool # for hybrid node, is it very bad triangle, udpate in updateGammaz!
    isBadTriangle::Bool # for hybrid node, is it very bad triangle, udpate in updateGammaz!
    inCycle::Int # = hybrid node if this node is part of a cycle created by such hybrid node, -1 if not part of cycle
    prev::Union{Nothing,ANode} # previous node in cycle, used in updateInCycle. set to "nothing" to begin with
    k::Int # num nodes in cycle, only stored in hybrid node, updated after node becomes part of network
           # default -1
    typeHyb::Int8 # type of hybridization (1,2,3,4, or 5), needed for quartet network only. default -1
    name::AbstractString
    clade::Vector{Int64} #Set to [1] initially to indicate that all start in main clade. outgroup added with addclades!
    # inner constructor: set hasHybEdge depending on edge
    Node() = new(-1,false,false,-1.,Edge[],false,false,false,false,false,false,-1,nothing,-1,-1,"",[1])
    Node(number::Int, leaf::Bool) = new(number,leaf,false,-1.,[],false,false,false,false,false,false,-1.,nothing,-1,-1,"",[1])
    Node(number::Int, leaf::Bool, hybrid::Bool) = new(number,leaf,hybrid,-1.,[],hybrid,false,false,false,false,false,-1.,nothing,-1,-1,"", [1])
    Node(number::Int, leaf::Bool, hybrid::Bool, edge::Array{Edge,1})=new(number,leaf,hybrid,-1.,edge,!all((e->!e.hybrid),edge),false,false,false,false,false,-1.,nothing,-1,-1,"", [1])
    Node(number::Int, leaf::Bool, hybrid::Bool,gammaz::Float64, edge::Array{Edge,1}) = new(number,leaf,hybrid,gammaz,edge,!all((e->!e.hybrid), edge),false,false,false,false,false,-1.,nothing,-1,-1,"", [1])
end

# partition type
mutable struct Partition
    cycle::Vector{Int} #hybrid node number for cycle (or cycles)
    edges::Vector{Edge} #edges in partition
end

abstract type Network end

# warning: no attempt to make sure the direction of edges matches with the root
# warning: no check if it is network or tree, node array can have no hybrids
# warning: nodes and edges need to be defined and linked before adding to a network
"""
    HybridNetwork

Subtype of abstract `Network` type.
Explicit network or tree with the following attributes:

- numTaxa
- numNodes (total number of nodes)
- numEdges
- numHybrids (number of hybrid nodes)
- edge (array of Edges)
- node (array of Nodes)
- root (index of root in vector 'node'. May be artificial, for printing and traversal purposes only.)
- hybrid (array of Nodes: those are are hybrid nodes)
- leaf (array of Nodes: those that are leaves)
- loglik (negative log pseudolik after estimation)
- isRooted (true or false)
"""
mutable struct HybridNetwork <: Network
    numTaxa::Int  # cannot set in constructor for congruence
    numNodes::Int
    numEdges::Int
    node::Array{Node,1}
    edge::Array{Edge,1}
    root::Int # node[root] is the root node, default 1
    names::Array{String,1} # translate table for taxon names --but also includes hybrid names...
    hybrid::Array{Node,1} # array of hybrid nodes in network
    numHybrids::Int # number of hybrid nodes
    cladewiseorder_nodeIndex::Vector{Int} # index in 'node' for "cladewise" preorder in main tree
    visited::Array{Bool,1} # reusable array of booleans
    edges_changed::Array{Edge,1} # reusable array of edges
    nodes_changed::Array{Node,1} # reusable array of nodes. used for preorder traversal
    leaf::Array{Node,1} # array of leaves
    ht::Vector{Float64} # vector of parameters to optimize
    numht::Vector{Int} # vector of number of the hybrid nodes and edges in ht e.g. [3,6,8,...], 2 hybrid nodes 3,6, and edge 8 is the 1st identifiable
    numBad::Int # number of bad diamond I hybrid nodes, set as 0
    hasVeryBadTriangle::Bool # true if the network has extremely/very bad triangles that should be ignored
    index::Vector{Int} #index in net.edge, net.node of elements in net.ht to make updating easy
    loglik::Float64 # value of the min -loglik after optBL
    blacklist::Vector{Int} # reusable array of integers, used in afterOptBL
    partition::Vector{Partition} # to choose edges from a partition only to avoid intersecting cycles
    cleaned::Bool # attribute to know if the network has been cleaned after readm default false
    isRooted::Bool # to know if network is rooted, e.g. after directEdges! (which updates isChild1 of each edge)
    # inner constructor
    function HybridNetwork(node::Array{Node,1},edge::Array{Edge,1})
        hybrid=Node[];
        leaf=Node[];
        for n in node
            if n.hybrid push!(hybrid,n); end
            if n.leaf   push!(leaf,  n); end
        end
        new(size(leaf,1),size(node,1),size(edge,1),node,edge,1,[],hybrid,size(hybrid,1), #numTaxa,...,numHybrids
            [],[],[],[],leaf,[],[], #cladewiseorder,...,numht
            0,false,[],0,[],[],false,false) #numBad...
    end
    function HybridNetwork(node::Array{Node,1},edge::Array{Edge,1},root::Int)
        hybrid=Node[];
        leaf=Node[];
        for n in node
            if n.hybrid push!(hybrid,n); end
            if n.leaf   push!(leaf,  n); end
        end
        new(size(leaf,1),size(node,1),size(edge,1),node,edge,root,[],hybrid,size(hybrid,1), #numTaxa,...,numHybrids
            [],[],[],[],leaf,[],[], #cladewiseorder,...,numht
            0,false,[],0,[],[],false,false) #numBad...
    end
    HybridNetwork() = new(0,0,0,[],[],0,[],[],0, # numTaxa ... numHybrid
                          [],[],[],[],[],[],[], # cladewiseorder...
                          0,false,[],0,[],[],false,false); # numBad ...
end

# type created from a HybridNetwork only to extract a given quartet
"""
    QuartetNetwork(net::HybridNetwork)

Subtype of `Network` abstract type.
need documentation!
"""
mutable struct QuartetNetwork <: Network
    numTaxa::Int
    numNodes::Int
    numEdges::Int
    node::Array{Node,1}
    edge::Array{Edge,1}
    hybrid::Array{Node,1} # array of hybrid nodes in network
    leaf::Array{Node,1} # array of leaves
    numHybrids::Int # number of hybrid nodes
    hasEdge::Array{Bool,1} # array of boolean with all the original identifiable edges of HybridNetwork and gammas (net.ht)
    quartetTaxon::Array{String,1} # the quartet taxa in the order it represents. Points to same array as its Quartet.taxon
    which::Int8 # 0 it tree quartet, 1 is equivalent to tree quartet and 2 if two minor CF different, default -1
    typeHyb::Array{Int8,1} #array with the type of hybridization of each hybrid node in the quartet
    t1::Float64 # length of internal edge, used when qnet.which=1, default = -1
    names::Array{String,1} # taxon and node names, same order as in network.node
    split::Array{Int8,1} # split that denotes to which side each leaf is from the split, i.e. [1,2,2,1] means that leaf1 and 4 are on the same side of the split, default -1,-1,-1,-1
    formula::Array{Int8,1} # array for qnet.which=1 that indicates if the expCf is major (1) or minor (2) at qnet.expCF[i] depending on qnet.formula[i], default -1,-1,-1
    expCF::Array{Float64,1} # three expected CF in order 12|34, 13|24, 14|23 (matching obsCF from qnet.quartet), default [0,0,0]
    indexht::Vector{Int} # index in net.ht for each edge in qnet.ht
    changed::Bool # true if the expCF would be changed with the current parameters in the optimization, to recalculate, default true
    index::Vector{Int} # index in qnet.edge (qnet.node for gammaz) of the members in qnet.indexht to know how to find quickly in qnet
    # inner constructor
    function QuartetNetwork(net::HybridNetwork)
        net2 = deepcopy(net); #fixit: maybe we dont need deepcopy of all, maybe only arrays
        new(net2.numTaxa,net2.numNodes,net2.numEdges,net2.node,net2.edge,net2.hybrid,net2.leaf,net2.numHybrids, [true for e in net2.edge],[],-1,[], -1.,net2.names,Int8[-1,-1,-1,-1],Int8[-1,-1,-1],[0,0,0],[],true,[])
        #new(sum([n.leaf?1:0 for n in net.node]),size(net.node,1),size(net.edge,1),copy(net.node),copy(net.edge),copy(net.hybrid),size(net.hybrid,1), [true for e in net2.edge],[],-1,[],-1.,net2.names,[-1,-1,-1,-1],[-1,-1,-1],[],true,[])
    end
    function QuartetNetwork(net::HybridNetwork,quartet::Array{String,1})
        net2 = deepcopy(net);
        new(net2.numTaxa,net2.numNodes,net2.numEdges,net2.node,net2.edge,net2.hybrid,net2.leaf,net2.numHybrids, [true for e in net2.edge],quartet,-1,[],-1.,net2.names,Int8[-1,-1,-1,-1],Int8[-1,-1,-1],[0,0,0],[],true,[])
    end
    QuartetNetwork() = new(0,0,0,[],[],[],[],0,[],[],-1,[],-1.0,[],[],[],[],[],true,[])
end

abstract type AQuartet end

"""
    Quartet

type that saves the information on a given 4-taxon subset. It contains the following attributes:

- number: integer
- taxon: vector of taxon names, like t1 t2 t3 t4
- obsCF: vector of observed CF, in order 12|34, 13|24, 14|23
- logPseudoLik
- ngenes: number of gene trees used to compute the observed CF; -1.0 if unknown
- qnet: [`QuartetNetwork`](@ref), which saves the expCF after snaq estimation to
  emphasize that the expCF depend on a specific network, not the data

see also: [`QuartetT`](@ref) for quartet with data of user-defined type `T`,
using a mapping between quartet indices and quartet taxa.
"""
mutable struct Quartet <: AQuartet
    number::Int
    taxon::Array{String,1} # taxa 1234. qnet.quartetTaxon points to the same array.
    obsCF::Array{Float64,1} # three observed CF in order 12|34, 13|24, 14|23
    qnet::QuartetNetwork # quartet network for the current network (want to keep as if private attribute)
    logPseudoLik::Float64 # log pseudolik value for the quartet. 0.0 by default
    ngenes::Float64 # number of gene trees used to compute the obsCV, default -1.; Float in case ngenes is average
    # inner constructor: to guarantee obsCF are only three and add up to 1
    function Quartet(number::Integer,t1::AbstractString,t2::AbstractString,t3::AbstractString,t4::AbstractString,obsCF::Array{Float64,1})
        size(obsCF,1) != 3 ? error("observed CF vector should have size 3, not $(size(obsCF,1))") : nothing
        0.99 < sum(obsCF) < 1.02 || @warn "observed CF should add up to 1, not $(sum(obsCF))"
        new(number,[t1,t2,t3,t4],obsCF,QuartetNetwork(),0.0,-1.0);
    end
    function Quartet(number::Integer,t1::Array{String,1},obsCF::Array{Float64,1})
        size(obsCF,1) != 3 ? error("observed CF vector should have size 3, not $(size(obsCF,1))") : nothing
        0.99< sum(obsCF) < 1.02 || @warn "observed CF should add up to 1, not $(sum(obsCF))"
        size(t1,1) != 4 ? error("array of taxa should have size 4, not $(size(t1,1))") : nothing
        0.0 <= obsCF[1] <= 1.0 || error("obsCF must be between (0,1), but it is $(obsCF[1]) for $(t1)")
        0.0 <= obsCF[2] <= 1.0 || error("obsCF must be between (0,1), but it is $(obsCF[2]) for $(t1)")
        0.0 <= obsCF[3] <= 1.0 || error("obsCF must be between (0,1), but it is $(obsCF[3]) for $(t1)")
        new(number,t1,obsCF,QuartetNetwork(),0.0,-1.0);
    end
    Quartet() = new(0,[],[],QuartetNetwork(),0.0,-1.0)
end

"""
QuartetT{T}

Generic type for 4-taxon sets. Fields:
- `number`: rank of the 4-taxon set
- `taxonnumber`: static vector of 4 integers, assumed to be distinct and sorted
- `data`: object of type `T`

For easier look-up, a unique mapping is used between the rank (`number`) of a
4-taxon set and its 4 taxa (see [`quartetrank`](@ref) and [`nchoose1234`](@ref)):

rank-1 = (t1-1) choose 1 + (t2-1) choose 2 + (t3-1) choose 3 + (t4-1) choose 4

# examples

```jldoctest
julia> nCk = PhyloNetworks.nchoose1234(5)
6×4 Array{Int64,2}:
 0   0   0  0
 1   0   0  0
 2   1   0  0
 3   3   1  0
 4   6   4  1
 5  10  10  5

julia> PhyloNetworks.QuartetT(1,3,4,6, [.92,.04,.04, 100], nCk)
4-taxon set number 8; taxon numbers: 1,3,4,6
data: [0.92, 0.04, 0.04, 100.0]
```
"""
struct QuartetT{T} <: AQuartet where T
    number::Int
    taxonnumber::StaticArrays.SVector{4,Int}
    data::T
end
function Base.show(io::IO, obj::QuartetT{T}) where T
    disp = "4-taxon set number $(obj.number); taxon numbers: "
    disp *= join(obj.taxonnumber,",")
    disp *= "\ndata: "
    print(io, disp)
    print(io, obj.data)
end
function QuartetT(tn1::Int,tn2::Int,tn3::Int,tn4::Int, data::T, nCk::Matrix, checksorted=true::Bool) where T
    if checksorted
        (tn1<tn2 && tn2<tn3 && tn3<tn4) || error("taxon numbers must be sorted")
    end
    QuartetT{T}(quartetrank(tn1,tn2,tn3,tn4,nCk), SVector(tn1,tn2,tn3,tn4), data)
end

"""
    quartetrank(t1,t2,t3,t4, nCk::Matrix)
    quartetrank([t1,t2,t3,t4], nCk)

Return the rank of a four-taxon set with taxon numbers `t1,t2,t3,t4`,
assuming that `ti`s are positive integers such that t1<t2, t2<t3 and t3<t4
(assumptions not checked!).
`nCk` should be a matrix of "n choose k" binomial coefficients:
see [`nchoose1234`](@ref).

# examples

```jldoctest
julia> nCk = PhyloNetworks.nchoose1234(5)
6×4 Array{Int64,2}:
 0   0   0  0
 1   0   0  0
 2   1   0  0
 3   3   1  0
 4   6   4  1
 5  10  10  5

julia> PhyloNetworks.quartetrank([1,2,3,4], nCk)
1

julia> PhyloNetworks.quartetrank([3,4,5,6], nCk)
15
```
"""
@inline function quartetrank(tnum::AbstractVector, nCk::Matrix)
    quartetrank(tnum..., nCk)
end
@inline function quartetrank(t1::Int, t2::Int, t3::Int, t4::Int, nCk::Matrix)
    # rank-1 = t1-1 choose 1 + t2-1 choose 2 + t3-1 choose 3 + t4-1 choose 4
    return nCk[t1,1] + nCk[t2,2] + nCk[t3,3] + nCk[t4,4] + 1
end

"""
    nchoose1234(nmax)

`nmax+1 x 4` matrix containing the binomial coefficient
"n choose k" in row `n+1` and column `k`. In other words,
`M[i,k]` gives "i-1 choose k". It is useful to store these
values and look them up to rank (a large number of) 4-taxon sets:
see [`quartetrank`](@ref).
"""
function nchoose1234(nmax::Int)
    # compute nC1, nC2, nC3, nC4 for n in [0, nmax]: used for ranking quartets
    M = Matrix{Int}(undef, nmax+1, 4)
    for i in 1:(nmax+1)
        M[i,1] = i-1 # n choose 1 = n. row i is for n=i-1
    end
    M[1,2:4] .= 0 # 0 choose 2,3,4 = 0
    for i in 2:(nmax+1)
        for k in 2:4 # to choose k items in 1..n: the largest could be n, else <= n-1
            M[i,k] = M[i-1,k-1] + M[i-1,k]
        end
    end
    return M
end

# Data on quartet concordance factors -------

"""
    DataCF

type that contains the following attributes:

- quartet (vector of Quartets)
- numQuartets
- tree (vector of trees: empty if a table of CF was input instead of list of trees)
- numTrees (-1 if a table CF was input instead of list of trees)
- repSpecies (taxon names that were repeated in table of CF or input gene trees: used inside snaq for multiple alleles case)

The list of Quartet may be accessed with the attribute .quartet.
If the input was a list of trees, the HybridNetwork's can be accessed with the attribute .tree.
For example, if the DataCF object is named d, d.quartet[1] will show the first quartet
and d.tree[1] will print the first input tree.
"""
mutable struct DataCF # fixit
    quartet::Array{Quartet,1} # array of quartets read from CF output table or list of quartets in file
    numQuartets::Integer # number of quartets
    tree::Vector{HybridNetwork} #array of input gene trees
    numTrees::Integer # number of gene trees
    repSpecies::Vector{String} #repeated species in the case of multiple alleles
    DataCF(quartet::Array{Quartet,1}) = new(quartet,length(quartet),[],-1,[])
    DataCF(quartet::Array{Quartet,1},trees::Vector{HybridNetwork}) = new(quartet,length(quartet),trees,length(trees),[])
    DataCF() = new([],0,[],-1,[])
end

# aux type for the updateBL function
mutable struct EdgeParts
    edgenum::Int
    part1::Vector{Node}
    part2::Vector{Node}
    part3::Vector{Node}
    part4::Vector{Node}
end

struct RootMismatch <: Exception
    msg::String
end
RootMismatch() = RootMismatch("")
Base.showerror(io::IO, e::RootMismatch) = print(io, "RootMismatch: ", e.msg);


text:
__precompile__()

struct Population
    law::Function
    fertilities::Array{Function, 1}
    initial::Float64
end

logistic(f::Number, p::Number) = f * p * (1-p)

mandelbrot(f::Number, p::Number) = f + p^2

ext: jl


text:
# some convenient type-piracy for scalars
AD.gradient(ad::AD.AbstractBackend, f, x::Real) = AD.derivative(ad, f, x)
AD.hessian(ad::AD.AbstractBackend, f, x::Real) = first.(AD.hessian(ad, f∘first, [x]))

# ForwardDiffBackend is built-in ot AutomaticDifferentiation
@init @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
    @eval const ForwardDiffBackend = AD.ForwardDiffBackend
end

# Zygote not yet in AutomaticDifferentiation, but remove this once it is
@init @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
    struct ZygoteBackend <: AD.AbstractBackend end
    AD.gradient(ad::ZygoteBackend, f, x) = Zygote.gradient(f, x)
    function AD.value_and_gradient(::ZygoteBackend, f, x)
        y, back = Zygote.pullback(f, x)
        y, back(Zygote.sensitivity(y))
    end
end


text:
@testset "Non-Trainable Potential Default Implmentation Unit Tests" begin
    p = MockNonTrainablePotential()

    @test get_parameters(p) == (;)
    @test set_parameters(p, (;)) === p
    @test serialize_parameters(p) == []
    @test deserialize_parameters(p, []) === p

    @test get_hyperparameters(p) == (;)
    @test set_hyperparameters(p, (;)) === p
    @test serialize_hyperparameters(p) == []
    @test deserialize_hyperparameters(p, []) === p
end


text:
module Kernels

using Distances
import LinearAlgebra

export SqExpKernel, metric

abstract type DistKernel end

struct SqExpKernel <: DistKernel
  alpha
  rho
  eps
end

SqExpKernel(alpha::Real, rho::Real) = SqExpKernel(alpha, rho, 1e-8)
metric(k::SqExpKernel) = SqEuclidean()

function (k::SqExpKernel)(Dsq::AbstractMatrix; jitter=true)
  out = k.alpha^2 * exp.(-Dsq / (2 * k.rho^2))
  if jitter
    out += LinearAlgebra.I * k.eps
  end
  return out
end

end

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
#= -------------------------------------------------------------------------
# @file problema1_con_disipador_parametric.jl
#
# @date 12/06/17 16:18:05
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
# imports
using DifferentialEquations
using PyPlot
#=------------------------------------------------------------------------------
                  sistem parametric equations
------------------------------------------------------------------------------=#
family = @ode_def heat_simulation begin
   dT_c = (-T_c / C_c)*( 1 / (R_ca) + 1 / (R_jc) + 1 / (R_cs))  + (T_j / (R_jc * C_c)) + (T_s / (R_cs * C_c))+ (T_a / (R_ca * C_c))
   dT_j = (T_c/ (R_jc * C_j)) - (T_j / (C_j * R_jc)) + (step_power(t) / C_j)
   dT_s = (T_c / (R_cs * C_s)) - (T_s / C_s) * (1 / (R_cs) + 1 / (R_sa)) + (T_a / (R_sa * C_s))
end R_jc=>1.0 C_c=>6.8 R_ca=>30.0 C_j=>0.5 T_a=>30.0 C_s=>39.6 R_cs=>0.4 R_sa=>4.5


text:
# ------------------------------------------------------------------------------
# enumerated types
# ------------------------------------------------------------------------------

module EvaluationMode

"""
Basis evaluation mode for operator inputs and outputs

# Modes:
- `interpolation`:     values interpolated to quadrature points
- `gradient`:          derivatives evaluated at quadrature points
- `quadratureweights`: quadrature weights

# Example:
```@jldoctest
EvaluationMode.EvalMode

# output
Enum LFAToolkit.EvaluationMode.EvalMode:
interpolation = 0
gradient = 1
quadratureweights = 2
```
"""
@enum EvalMode interpolation gradient quadratureweights

end # submodule

module MultigridType

"""
Multigrid types

# Types:
- `pmultigrid`: p-multigrid
- `hmultigrid`: h-multigrid

# Example:
```@jldoctest
LFAToolkit.MultigridType.MgridType

# output
Enum LFAToolkit.MultigridType.MgridType:
pmultigrid = 0
hmultigrid = 1
```
"""
@enum MgridType hmultigrid pmultigrid

end # submodule

module BDDCInjectionType

"""
BDDC injection types

# Types:
- `scaled`:   scaled injection
- `harmonic`: discrete harmonic extension

# Example
```@jldoctest
LFAToolkit.BDDCInjectionType.BDDCInjectType

# output
Enum LFAToolkit.BDDCInjectionType.BDDCInjectType:
scaled = 0
harmonic = 1
```
"""
@enum BDDCInjectType scaled harmonic

end # submodule

# ------------------------------------------------------------------------------


text:
const K_PARAM = """    k::Integer
Number of neighbors (must be greater than 0)."""

const KNN_PARAMS = """    metric::Metric
This is one of the Metric types defined in the Distances.jl package. It is possible to define your own metrics by
creating new types that are subtypes of Metric.

    algorithm::Symbol
One of `(:kdtree, :balltree)`. In a `kdtree`, points are recursively split into groups using hyper-planes.
Therefore a KDTree only works with axis aligned metrics which are: Euclidean, Chebyshev, Minkowski and Cityblock.
A *brutetree* linearly searches all points in a brute force fashion and works with any Metric. A *balltree*
recursively splits points into groups bounded by hyper-spheres and works with any Metric.

    static::Union{Bool, Symbol}
One of `(true, false, :auto)`. Whether the input data for fitting and transform should be statically or dynamically
allocated. If `true`, the data is statically allocated. If `false`, the data is dynamically allocated. If `:auto`,
the data is dynamically allocated if the product of all dimensions except the last is greater than 100.

    leafsize::Int
Determines at what number of points to stop splitting the tree further. There is a trade-off between traversing the
tree and having to evaluate the metric function for increasing number of points.

    reorder::Bool
While building the tree this will put points close in distance close in memory since this helps with cache locality.
In this case, a copy of the original data will be made so that the original data is left unmodified. This can have a
significant impact on performance and is by default set to true.

    parallel::Bool
Parallelize `score` and `predict` using all threads available. The number of threads can be set with the
`JULIA_NUM_THREADS` environment variable. Note: `fit` is not parallel."""

# for some reason providing this as a function is massively slowing down compilation
macro tree(detector, X)
    esc(quote
        $detector.algorithm === :kdtree ?
            NN.KDTree($X, $detector.metric; $detector.leafsize, $detector.reorder) :
            NN.BallTree($X, $detector.metric; $detector.leafsize, $detector.reorder)
    end)
end

function knn_parallel(tree::NN.NNTree, X::AbstractArray, k::Int,
    sort::Bool = false)::Tuple{AbstractVector, AbstractVector}
    # pre-allocate the result arrays (as in NearestNeighbors.jl)
    samples = size(X, 2)
    dists = [Vector{NN.get_T(eltype(X))}(undef, k) for _ in 1:samples]
    idxs = [Vector{Int}(undef, k) for _ in 1:samples]

    # get number of threads
    nThreads = Threads.nthreads()
    # partition the input array equally
    partition_size = samples ÷ nThreads + 1
    partitions = Iterators.partition(axes(X, 2), partition_size)
    Threads.@threads for idx = collect(partitions)
        @inbounds idxs[idx], dists[idx] = NN.knn(tree, view(X, :, idx), k, sort)
    end
    idxs, dists
end

function dnn_parallel(tree::NN.NNTree, X::AbstractArray, d::Real, sort::Bool = false)::AbstractVector
    # pre-allocate the result arrays (as in NearestNeighbors.jl)
    samples = size(X, 2)
    scores = Vector{Float64}(undef, samples)

    # get number of threads
    nThreads = Threads.nthreads()
    # partition the input array equally
    partition_size = samples ÷ nThreads + 1
    partitions = Iterators.partition(axes(X, 2), partition_size)
    Threads.@threads for idx = collect(partitions)
        @inbounds scores[idx] = dnn(NN.inrange(tree, view(X, :, idx), d, sort))
    end
    scores
end

# Calculate the k-nearest neighbors with ignoring the own point in the tree.
function knn_others(tree::NN.NNTree, X::AbstractArray, k::Integer)::Tuple{AbstractVector, AbstractVector}
    idxs, dists = NN.knn(tree, X, k + 1, true) # we ignore the distance to the 'self' point, important to sort!
    ignore_self = vecvec -> map(vec -> vec[2:end], vecvec)
    ignore_self(idxs), ignore_self(dists)
end

# The NN package automatically converts matrices to a vector of points (static vectors) for improved performance
# this results in very bad performance for high-dimensional matrices (e.g. d > 100).
dynamic_view(X::Data) = [NN.SizedVector{length(v)}(v) for v in eachslice(X; dims=ndims(X))]
auto_view(X::Data) = prod(size(X)[1:end-1]) > 100 ? dynamic_view(X) : X
prepare_data(X::Data, static::Union{Bool, Symbol}) = static === :auto ? auto_view(X) :
                                                     static === false ? dynamic_view(X) : X


text:
"""
std_w

Weighted standard deviation.
Negative weights are not allowed.
"""
function std_w(xM :: Array{T}, wtM :: Array{T}) where T <: AbstractFloat
	xV = vec(xM)
	wtV = AnalyticWeights(vec(wtM))

	## Input check
	if any(wtV .< 0.0)
	  error("Invalid weights")
	end

	# wtSum = sum(wtV);
	# xMean = sum(xV .* wtV) / wtSum;
	# xStd  = sqrt(sum(((xV .- xMean) .^ 2) .* wtV)  /  wtSum);

	xMean, xStd = mean_and_std(xV, wtV, corrected = false);
	# xStd = std(xV, wtV);
	return xStd, xMean
end


text:
"""
    TrivialCallback()

A callback that does nothing. This can be useful to disable some callbacks
easily via [`trixi_include`](@ref).
"""
function TrivialCallback()
  DiscreteCallback(trivial_callback, trivial_callback,
                   save_positions=(false,false))
end

trivial_callback(u, t, integrator) = false
trivial_callback(integrator) = u_modified!(integrator, false)

function Base.show(io::IO, cb::DiscreteCallback{Condition,Affect!}) where {Condition, Affect!<:typeof(trivial_callback)}
  print(io, "TrivialCallback()")
end

# This allows to set `summary_callback = TrivialCallback()` in elixirs to suppress
# output, e.g. in `convergence_test`.
function (cb::DiscreteCallback{Condition,Affect!})(io::IO=stdout) where {Condition, Affect!<:typeof(trivial_callback)}
  return nothing
end


text:
"""
    evolve!(curr::Field, prev::Field, a, dt)

Calculate a new temperature field curr based on the previous
field prev. a is the diffusion constant and dt is the largest
stable time step.
"""
function evolve!(curr::Field, prev::Field, a, dt)
    for j = 2:curr.ny+1
        for i = 2:curr.nx+1
            xderiv = (prev.data[i-1, j] - 2.0 * prev.data[i, j] + prev.data[i+1, j]) / curr.dx^2
            yderiv = (prev.data[i, j-1] - 2.0 * prev.data[i, j] + prev.data[i, j+1]) / curr.dy^2
            curr.data[i, j] = prev.data[i, j] + a * dt * (xderiv + yderiv)
        end
    end
end

"""
    swap_fields!(curr::Field, prev::Field)

Swap the data of two fields curr and prev.
"""
function swap_fields!(curr::Field, prev::Field)
    tmp = curr.data
    curr.data = prev.data
    prev.data = tmp
end

"""
    average_temperature(f::Field)

Calculate average temperature of a temperature field.
"""
average_temperature(f::Field) = sum(f.data[2:f.nx+1, 2:f.ny+1]) / (f.nx * f.ny)

"""
    simulate!(current, previous, nsteps)

Run the heat equation solver on fields curr and prev for nsteps.
"""
function simulate!(curr::Field, prev::Field, nsteps)

    println("Initial average temperature: $(average_temperature(curr))")

    # Diffusion constant
    a = 0.5
    # Largest stable time step
    dt = curr.dx^2 * curr.dy^2 / (2.0 * a * (curr.dx^2 + curr.dy^2))

    # display a nice progress bar
    p = Progress(nsteps)

    for i = 1:nsteps
        # calculate new state based on previous state
        evolve!(curr, prev, a, dt)

        # swap current and previous fields
        swap_fields!(curr, prev)

        # increment the progress bar
        next!(p)
    end

    # print final average temperature
    println("Final average temperature: $(average_temperature(curr))")
end


text:
function constant_features_edge(s::Integer, d::Integer, mapf::MAPF)
    g = mapf.graph
    outdeg_s = outdegree(g, s)
    indeg_s = indegree(g, s)
    outdeg_d = outdegree(g, d)
    indeg_d = indegree(g, d)
    rev_edge_exists = has_edge(g, d, s)
    return Float64[outdeg_s, outdeg_d]
end

function solution_features_edge(s::Integer, d::Integer, solution::Solution, mapf::MAPF)
    g = mapf.graph
    paths_crossing_s = 0
    paths_crossing_d = 0
    paths_crossing_e = 0
    for path in solution
        for ((t1, v1), (t2, v2)) in zip(path[1:(end - 1)], path[2:end])
            paths_crossing_s += (v1 == s) + (v2 == s)
            paths_crossing_d += (v1 == d) + (v2 == d)
            paths_crossing_e += (v1 == s && v2 == d)
        end
    end
    return Float64[paths_crossing_s, paths_crossing_d]
end

function all_features_edge(s::Integer, d::Integer, solution::Solution, mapf::MAPF)
    return vcat(
        constant_features_edge(s, d, mapf), solution_features_edge(s, d, solution, mapf)
    )
end

function edges_embedding(mapf::MAPF)
    solution = independent_dijkstra(mapf)
    x = reduce(
        hcat,
        all_features_edge(src(ed), dst(ed), solution, mapf) for ed in edges(mapf.graph)
    )
    s = std(x; dims=2)
    s[isapprox.(s, 0.0)] .= 1  # columns with zero variance
    x = (x .- mean(x; dims=2)) ./ s
    return x
end


text:
using DashBootstrapComponents, DashHtmlComponents

buttons = html_div([
    dbc_button("Primary", outline=true, color="primary", className="mr-1"),
    dbc_button("Secondary", outline=true, color="secondary", className="mr-1"),
    dbc_button("Success", outline=true, color="success", className="mr-1"),
    dbc_button("Warning", outline=true, color="warning", className="mr-1"),
    dbc_button("Danger", outline=true, color="danger", className="mr-1"),
    dbc_button("Info", outline=true, color="info", className="mr-1"),
    dbc_button("Light", outline=true, color="light", className="mr-1"),
    dbc_button("Dark", outline=true, color="dark"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
function _conic_branch!(t::ExtensionType, m::GlobalOptimizer{N,T}) where {N,T}
    return false
end


text:
"""
    ArrayNode{A <: AbstractArray, C} <: AbstractMillNode

Data node for storing array-like data of type `A` and metadata of type `C`. The convention is that samples are stored along the last axis, e.g. in columns of a matrix.

See also: [`AbstractMillNode`](@ref), [`ArrayModel`](@ref).
"""
struct ArrayNode{A <: AbstractArray, C} <: AbstractMillNode
    data::A
    metadata::C
end

"""
    ArrayNode(d::AbstractArray, m=nothing)

Construct a new [`ArrayNode`](@ref) with data `d` and metadata `m`.

# Examples
```jldoctest
julia> a = ArrayNode([1 2; 3 4; 5 6])
3×2 ArrayNode{Matrix{Int64}, Nothing}:
 1  2
 3  4
 5  6
```

See also: [`AbstractMillNode`](@ref), [`ArrayModel`](@ref).
"""
ArrayNode(d::AbstractArray) = ArrayNode(d, nothing)

Flux.@functor ArrayNode

mapdata(f, x::ArrayNode) = ArrayNode(mapdata(f, x.data), x.metadata)

dropmeta(x::ArrayNode) = ArrayNode(x.data)

Base.ndims(x::ArrayNode) = Colon()
Base.size(x::ArrayNode) = size(x.data)

StatsBase.nobs(a::ArrayNode) = size(a.data, 2)
StatsBase.nobs(a::ArrayNode, ::Type{ObsDim.Last}) = nobs(a)

function reduce(::typeof(catobs), as::Vector{<:ArrayNode})
    ArrayNode(reduce(catobs, data.(as)), reduce(catobs, metadata.(as)))
end

_cat_meta(f, m::Vector{Nothing}) = nothing
_cat_meta(f, m) = reduce(f, m)

Base.vcat(as::ArrayNode...) = reduce(vcat, collect(as))
function Base.reduce(::typeof(vcat), as::Vector{<:ArrayNode})
    ArrayNode(reduce(vcat, data.(as)), _cat_meta(vcat, metadata.(as)))
end

Base.hcat(as::ArrayNode...) = reduce(hcat, collect(as))
function Base.reduce(::typeof(hcat), as::Vector{<:ArrayNode})
    ArrayNode(reduce(hcat, data.(as)), _cat_meta(hcat, metadata.(as)))
end

Base.getindex(x::ArrayNode, i::VecOrRange{<:Int}) = ArrayNode(subset(x.data, i), subset(x.metadata, i))

Base.hash(e::ArrayNode, h::UInt) = hash((e.data, e.metadata), h)
(e1::ArrayNode == e2::ArrayNode) = isequal(e1.data == e2.data, true) && e1.metadata == e2.metadata
Base.isequal(e1::ArrayNode, e2::ArrayNode) = isequal(e1.data, e2.data) && isequal(e1.metadata, e2.metadata)

function Base.show(io::IO, ::MIME"text/plain", @nospecialize(n::ArrayNode))
    print(io, join(size(n.data), "×"), " ", summary(n))
    if !isempty(n.data)
        print(io, ":\n")
        Base.print_array(IOContext(io, :typeinfo => eltype(n.data)), n.data)
    end
end

function _show_data(io, n::ArrayNode{T}) where T <: AbstractArray
    print(io, "(")
    if ndims(n.data) == 1
        print(io, nameof(T), " of length ", length(n.data))
    else
        print(io, join(size(n.data), "×"), " ", nameof(T))
    end
    print(io, " with ", eltype(n.data), " elements)")
end


text:
using Test
using Z3

ctx = Context()
x = real_const(ctx, "x")
y = real_const(ctx, "y")

s = Solver(ctx, "QF_NRA")
add(s, x == y^2)
add(s, x > 1)

res = check(s)
@test res == Z3.sat


text:
import PrettyTables: pretty_table

mutable struct Wave
    ω::Real
    β::Real
    n1::Number
    n2::Number
    n′::Number
end

mutable struct Setup
    Ω::Wave
    Ω1::Wave
    Ω2::Wave
end

"""
    Setup(ω, β, n1, n2)

Create an object representing the setup.

# Arguments
* `ω`: `Tuple{Real,Real}` containing \$\\omega_1\$ and \$\\omega_2\$ so that \$\\omega = \\omega_1 + \\omega_2\$.
* `β`: `Tuple{Real,Real}` containing the incident angles corresponing to \$\\omega_1\$ and \$\\omega_2\$.
* `n1`: `Tuple{Number,Number,Number}` containing the refractive indices of medium 1 corresponding to \$\\omega\$, \$\\omega_1\$ and \$\\omega_2\$.
* `n2`: `Tuple{Number,Number,Number}` containing the refractive indices of medium 2 corresponding to \$\\omega\$, \$\\omega_1\$ and \$\\omega_2\$.
"""
function Setup(
        (ω1, ω2)::Tuple{Real,Real},
        (β1, β2)::Tuple{Real,Real},
        (n1, n11, n12)::Tuple{Number,Number,Number},
        (n2, n21, n22)::Tuple{Number,Number,Number},
    )
    ω = ω1 + ω2
    β = sfangle(ω, ω1, ω2, β1, β2)
    n′, n′1, n′2 = ninterface.((n2, n21, n22))

    wave  = Wave(ω, β, n1, n2, n′)
    wave1 = Wave(ω1, β1, n11, n21, n′1)
    wave2 = Wave(ω2, β2, n12, n22, n′2)

    Setup(wave, wave1, wave2)
end

function Base.show(io::IO, s::Setup)
    ω = s.Ω.ω
    ω1= s.Ω1.ω
    ω2= s.Ω2.ω
    λ = freq2wl(ω)
    λ1= freq2wl(ω1)
    λ2= freq2wl(ω2)

    data = [
        "Ω"  λ  ω  rad2deg(s.Ω.β)  s.Ω.n1  s.Ω.n2  s.Ω.n′;
        "Ω1" λ1 ω1 rad2deg(s.Ω1.β) s.Ω1.n1 s.Ω1.n2 s.Ω1.n′;
        "Ω2" λ2 ω2 rad2deg(s.Ω2.β) s.Ω2.n1 s.Ω2.n2 s.Ω2.n′;
    ]
    formatter = (v,i,j) -> typeof(v) <: Number ? round(v, digits=3) : v
    header = ["", "λ (nm)", "ω (s⁻¹)", "β (deg)", "n₁", "n₂", "n′"]

    println("SFGAnalysis.Setup")
    pretty_table(data; formatters=formatter, header=header)
end

# In order to not broadcast Setup and treat it as a scalar
Broadcast.broadcastable(s::Setup) = (s,)


text:
module FLuxSymbolEquationsV2
# excitvar à automatiser OK
# affectation des valeurs pour les composants
# afectation des variables pour plot
# variables d'intégration automatique OK
# just for test
export SolveSymbolic,NumericSolve,PlotSol

using Plots
using DifferentialEquations
using SymPy
mutable struct Schematics
	Template::Array{String}
end

mutable struct Component

	name::String
	loopNb::Int64
	formula::String
	position::String
	current::String
	special::Bool
	ctype::String
	value::Float64
	rect::Bool
	commons::Int64
	voltage::String
	Component()=new()
end

mutable struct LoopC
	CompInLoop::Array{Component}
	CommonComp::Array{Component}
	LoopNumber::Int64
	LoopC()=new()
end

mutable struct CommonComp
	Common::Component
	Loops::Array{Int64}
end

myInputCircuit=["!VUd 110.,Rp 1e-4 ,Lr 32.4e-6,T1P 1.82e-3,MT1 0.000823308666297155,Cr 39e-9","MT1 ,T1S 380e-6 ,R2 1e-4 ,!rect ,Co 10e-6","Co ,Rl 2.4"]

#export greet, myStrgToSolve, myGraph

function formul(typ,nam,cur)
	formul=""
	if typ=="R"
		formul="V"*nam*"="*nam*"*"*cur
	elseif typ=="L"
		formul="V"*nam*"="*nam*"*d"*cur
	elseif typ=="C"
		formul=cur*"="*nam*"*dV"*nam
	elseif typ=="T"
		formul="V"*nam*"="*nam*"*d"*cur
	elseif typ=="M"
		formul="V"*nam*"=-"*nam*"*d"*cur
	elseif typ=="m"
		formul="V"*nam*"=-"*"M"*nam[2:end]*"*d"*cur
	end
	return formul
end
# flatten formulas
function flattenEqs(LLCExample,CommonCmpt,myLoop1)
	flattenEq=""
	# process voltage sum in loop equal zero
	sumvolt=String[]
	for i in 1:length(LLCExample.Template ) push!(sumvolt,"") end
	for i in 1:length(myLoop1)
		sumvolt[myLoop1[i].loopNb]*=myLoop1[i].voltage*"+"
	end
	for j in 1:length(LLCExample.Template)
		sumvolt[j]=sumvolt[j][1:end-1]
	end
	for j in 1:4:length(CommonCmpt)
		if CommonCmpt[j][1:1]!="M"
			sumvolt[CommonCmpt[j+3]]=replace(sumvolt[CommonCmpt[j+3]],"V"*CommonCmpt[j]=>"(-V"*CommonCmpt[j]*")")
		end
	end
	for i in 1:length(sumvolt)
		sumvolt[i]=replace(sumvolt[i],"++"=>"+")
	end
	# delete extra formula
	for i in length(myLoop1):-1:1
		for j in 1:4:length(CommonCmpt)
			if (myLoop1[i].name ==CommonCmpt[j]) && myLoop1[i].current=="I"*string(CommonCmpt[j+3])  && myLoop1[i].ctype!="M"
				splice!(myLoop1, i)
			end
		end
	end

	# build list of equation
	for i in 1:length(myLoop1)
		flattenEq*=replace(myLoop1[i].formula,"="=>"-(")*"),"
	end
	flattenEq=flattenEq[3:end]
	for i in 1:length(sumvolt)
		sumvolt[i]=replace(sumvolt[i],"++"=>"+")
		flattenEq*=sumvolt[i]*","
	end

	flattenEq=replace(flattenEq,",),"=>",")

	return flattenEq[1:end-1],sumvolt
end

function fillCmpt(LLCExample)
	rectOn=false
	myLoop=[]
	CommonCmpt=[]
	CommonLoop=[]
	test=[]
    excitVar0=[]
	for j in 1:length(LLCExample.Template)  # numero de boucle
		parts=split(LLCExample.Template[j],",")
		for i in 1:length(parts)
			tCmpt=Component()
			temp0=split(parts[i]," ")
			if temp0[1][1:1]=="!"
				tCmpt.special=true
				tCmpt.name=temp0[1][2:end]
				if tCmpt.name=="rect"
					rectOn=true
					tCmpt.voltage=""
				else
					tCmpt.voltage=tCmpt.name
					push!(excitVar0,tCmpt.name)
				end

			else
				tCmpt.special=false
				tCmpt.name=temp0[1][1:end]

				if temp0[2]==""  !!!!!!!!!!!!!!!
					tCmpt.commons=j
					push!(CommonCmpt,tCmpt.name)
					push!(CommonLoop,j)
					#tCmpt.value=1e-139 # dummy value
				end

				tCmpt.voltage="V"*tCmpt.name

			end
			tCmpt.ctype=tCmpt.name[1:1]
			tCmpt.loopNb=j
			tCmpt.current="I"*string(tCmpt.loopNb)
			tCmpt.formula=formul(tCmpt.ctype,tCmpt.name,tCmpt.current)
			tCmpt.rect=rectOn
			#if tCmpt.value!=false
			try
				tCmpt.value=parse(Float64, temp0[2])
			catch
				for k in 1:length(myLoop)
					if myLoop[k].name==tCmpt.name
						tCmpt.value=myLoop[k].value
						break
					end
				end
			end
			#end
			push!(myLoop,tCmpt)
		end
	end
	for i in 1:length(myLoop)
		if myLoop[i].name in CommonCmpt
			push!(test,myLoop[i].name)
			push!(test,myLoop[i].loopNb)
		end
	end
	return myLoop,test,excitVar0
end
function proct(myLoop,CommonCmpt)
	cur=[]
	for i in 1:length(myLoop)
		if myLoop[i].ctype=="M"
			push!(cur,myLoop[i].current)
		end
	end
	fstok=0
	for i in 1:length(myLoop)
		if myLoop[i].ctype=="M"	&& fstok==0
			myLoop[i].current=pop!(cur)
			#myLoop[i].formula=myLoop[i].formula[1:end-2]*myLoop[i].current
			myLoop[i].formula=myLoop[i].formula[1:4]*"P"*myLoop[i].formula[5:end-2]*myLoop[i].current

			myLoop[i].voltage=myLoop[i].voltage*"P"
			fstok=1
		elseif myLoop[i].ctype=="M"	&& fstok==1
			myLoop[i].current=pop!(cur)
			#myLoop[i].formula=myLoop[i].formula[1:end-2]*myLoop[i].current
			myLoop[i].formula=myLoop[i].formula[1:4]*"S"*myLoop[i].formula[5:end-2]*myLoop[i].current
			myLoop[i].voltage=myLoop[i].voltage*"S"
			fstok=0
		end
	end

	for i in 1:2:length(CommonCmpt)   # courants impliqués  cur= [2, 3]
		if CommonCmpt[i][1:1]!="M"
			push!(cur,CommonCmpt[i+1])
		end
	end
	for i in 1:4:length(CommonCmpt)
		if CommonCmpt[i][1:1]!="M"
			for j in 1:length(myLoop)    # process for common componzents
				if myLoop[j].name==CommonCmpt[i]   #replace I2  with (I2-I3)
					myLoop[j].formula=replace(myLoop[j].formula,"I"*string(cur[1])=>"(I"*string(cur[1])*"-I"*string(cur[2])*")")
				end
			end
		end
	end
	return  myLoop
end
function getmyVars(myLoop1,LLCExample)
	mySymbvars=[]  # all symbolic var
	myvars=[]   # var to be integrated by ODE
	myvars1=""
	for i in 1:length(myLoop1)
		if myLoop1[i].ctype in ["L","M","T"]
			push!(myvars,"d"*myLoop1[i].current)
			push!(mySymbvars,myLoop1[i].current)
		elseif myLoop1[i].ctype=="C"
			push!(myvars,"d"*myLoop1[i].voltage)
			push!(mySymbvars,myLoop1[i].voltage)

		end
	end
	for i in 1:length(myLoop1)
		if myLoop1[i].name !="" && myLoop1[i].name !="rect"
		#push!(myvars, myLoop1[i].current)
			if myLoop1[i].special==false && ("d"*myLoop1[i].voltage in myvars) == false
				push!(myvars,myLoop1[i].voltage)
			else
				push!(mySymbvars,myLoop1[i].voltage)
			end
		push!(mySymbvars,myLoop1[i].name)
		end
	end
	myvars=unique(myvars)
    mySymbvars=unique(mySymbvars)
	#il faut rajouter I3 car pas dI3
	for i in 1:length(LLCExample.Template)
		if ("dI"*string(i) in myvars) == false
			push!(myvars,"I"*string(i))
		end
	end
	return myvars,mySymbvars
end
function flatVars(myVar)
	myvars1=""
   for i in 1:length(myVar[1])
		myvars1*=myVar[1][i]*","
	end
	return myvars1[1:end-1]
end
# create rectComp="I2,Co,Rl"
function rectdata(myLoop1)
	rectComp=""
	for i in 1:length(myLoop1)
		if myLoop1[i].name=="rect"
			rectComp*=myLoop1[i].current*","
		elseif myLoop1[i].rect==true
			rectComp*=myLoop1[i].name*","
		end
	end
		return rectComp[1:end-1]
end

function rect(rectComp,mydifEq)
	compsRect=split(rectComp,",")
	cour=""
	if compsRect!=[] cour=compsRect[1] end
	my=[]
	for i in 1:length(compsRect)
		if compsRect[i][1:1]=="C"
			for j in 1:length(mydifEq)
				if mydifEq[j][1:2+length(compsRect[i])]=="dV"*compsRect[i]
					push!(my,replace(mydifEq[j],cour=>"abs("*cour*")"))
				else #Uco*(-sign(Is))
					push!(my,replace(mydifEq[j],"V"*compsRect[i]=>"V"*compsRect[i]*"*(sign("*cour*"))"))
				end
			end
		end
	end

	return my,cour
end

function collectVar(rt)
	myst=""
	myVrArray=[]
	for i in 1:length(rt)
	myst=myst*rt[i][2:findfirst(isequal('='), rt[i])-1]*","
	push!(myVrArray,rt[i][2:findfirst(isequal('='), rt[i])-1])
	end
	return myst[1:end-1],myVrArray
end

function flatEqs(sympyEq,myVar)
	eqstoSolv=""
	myflatVar=""
	for ii in 1:length(sympyEq)-1
		eqstoSolv=eqstoSolv*sympyEq[ii]*","
	end
	eqstoSolv=eqstoSolv*sympyEq[length(sympyEq)]
	for ii in 1:length(myVar)-1
		myflatVar=myflatVar*myVar[1][ii]*","
	end
	myflatVar=myflatVar*myVar[1][length(myVar)]
	return eqstoSolv,myflatVar
end

# function giving equation from circuit
function myEquationsSet(myInputCircuit)
	mys=""
	excitVar=""  #
	LLCExample=Schematics(myInputCircuit)
	myLoop,CommonCmpt,excitVar=fillCmpt(LLCExample)
	myLoop1=proct(myLoop,CommonCmpt)
	flattenEqs1,sumvolt=flattenEqs(LLCExample,CommonCmpt,myLoop1)
	myVar=getmyVars(myLoop1,LLCExample)
	symbdat=""
	#Création des noms symboliques)
	for j in 1:length(myVar)
		for i in 1:length(myVar[j])
			try
				symbdat=replace(replace(replace(string(myVar[j][i]),"\"" =>""),"Any["=>""),"]"=>"")
				eval(Meta.parse(symbdat*"=symbols(\""*symbdat*"\")"))
			catch
			end
		end
	end

	Az=linsolve(eval(Meta.parse(flattenEqs1)),eval(Meta.parse(flatVars(myVar))))
	mysol=split(replace(string.(Az),"FiniteSet(("=>"")[1:end-2],",")
	VartoSolve=myVar[1]
	mydVoudIindx=findall.("d",VartoSolve)
	difEqIndex=findall(mydVoudIindx -> mydVoudIindx!=[],mydVoudIindx)
    mydifEq=VartoSolve[difEqIndex].*"=".*mysol[difEqIndex]
	rectComp=rectdata(myLoop1)
	compsRect=split(rectComp,",")
	cour=rect(rectComp,mydifEq)[2]
	if compsRect!=[] cour=compsRect[1] end
    myst,myVrArray=collectVar(rect(rectComp,mydifEq)[1])
	myStrgToSolve=""
	myStrgToSolve1=""
	rect1=[]
	for i in 1:length(rect(rectComp,mydifEq)[1])
		push!(rect1,"du["*string(i)*"] = "*rect(rectComp,mydifEq)[1][i])
	end
	excitVar=fillCmpt(LLCExample)[3]
	#replace.(rect1,excitVar=>excitVar*"*p[1]")
	myStrgToSolve=flatEqs(replace.(rect1,excitVar[1]=>excitVar[1]*"*p[1]"),myVar)[1]*"; end"
	myStrgToSolve1="function eqDiffToSolve!(du,u,p,t) "* myst *" = u; " *replace(myStrgToSolve,","=>"; ") #!!!!!!
	eval(Meta.parse(myStrgToSolve1))
	freq=1e5  #OK
	u₀=[0,0,0,0]
	p=[-1.0,1.0]
	tspan = (0.0,600.e-6)
	nbPeriod=60
	dosetimes = collect(1:nbPeriod)/freq/2
	#eval(Meta.parse(myStrgToSolve1))
	for i in 1:length(myLoop)
		#if typeof(myLoop[i].value)!=1e-139
			try eval(Meta.parse(myLoop[i].name*"="*string(myLoop[i].value)))
			catch
			end
		#end
	end
	prob = ODEProblem(eqDiffToSolve! ,u₀,tspan,p)
	condition(u,t,integrator) = t ∈ dosetimes
	affect!(integrator) = integrator.p[1] = -integrator.p[1]
	cb = DiscreteCallback(condition,affect!)
	sol= DifferentialEquations.solve(prob ,Tsit5(),callback=cb,tstops=dosetimes)
    myGraph=plot(sol,linewidth=2,xaxis="t",label=permutedims(myVrArray),layout=(4,1))
	return myGraph
    #return myStrgToSolve1
	#return excitVar,myStrgToSolve
	#return rect1
    #return myst,myVrArray
	#return cour
	#return VartoSolve.*"=".*mysol
	#return VartoSolve[difEqIndex].*"=".*mysol[difEqIndex]
end

function SolveSymbolic(myInputCircuit)
	mys=""
	excitVar=""  #
	LLCExample=Schematics(myInputCircuit)
	myLoop,CommonCmpt,excitVar=fillCmpt(LLCExample)
	myLoop1=proct(myLoop,CommonCmpt)
	flattenEqs1,sumvolt=flattenEqs(LLCExample,CommonCmpt,myLoop1)
	myVar=getmyVars(myLoop1,LLCExample)
	symbdat=""
	#Création des noms symboliques)
	for j in 1:length(myVar)
		for i in 1:length(myVar[j])
			try
				symbdat=replace(replace(replace(string(myVar[j][i]),"\"" =>""),"Any["=>""),"]"=>"")
				eval(Meta.parse(symbdat*"=symbols(\""*symbdat*"\")"))
			catch
			end
		end
	end

	Az=linsolve(eval(Meta.parse(flattenEqs1)),eval(Meta.parse(flatVars(myVar))))
	mysol=split(replace(string.(Az),"FiniteSet(("=>"")[1:end-2],",")
	VartoSolve=myVar[1]
	mydVoudIindx=findall.("d",VartoSolve)
	difEqIndex=findall(mydVoudIindx -> mydVoudIindx!=[],mydVoudIindx)
    mydifEq=VartoSolve[difEqIndex].*"=".*mysol[difEqIndex]
	rectComp=rectdata(myLoop1)
	compsRect=split(rectComp,",")
	cour=rect(rectComp,mydifEq)[2]
	if compsRect!=[] cour=compsRect[1] end
    myst,myVrArray=collectVar(rect(rectComp,mydifEq)[1])
	myStrgToSolve=""
	myStrgToSolve1=""
	rect1=[]
	for i in 1:length(rect(rectComp,mydifEq)[1])
		push!(rect1,"du["*string(i)*"] = "*rect(rectComp,mydifEq)[1][i])
	end
	excitVar=fillCmpt(LLCExample)[3]
	#replace.(rect1,excitVar=>excitVar*"*p[1]")
	myStrgToSolve=flatEqs(replace.(rect1,excitVar[1]=>excitVar[1]*"*p[1]"),myVar)[1]*"; end"
	myStrgToSolve1="function eqDiffToSolve!(du,u,p,t) "* myst *" = u; " *replace(myStrgToSolve,","=>"; ") #!!!!!!
	eval(Meta.parse(myStrgToSolve1))
	return myStrgToSolve1,myLoop,myVrArray
end

function NumericSolve(myStrgToSolve1,myLoop)
	freq=1e5  #OK
	u₀=[0,0,0,0]
	p=[-1.0,1.0]
	tspan = (0.0,600.e-6)
	nbPeriod=60
	dosetimes = collect(1:nbPeriod)/freq/2
	#eval(Meta.parse(myStrgToSolve1))
	for i in 1:length(myLoop)
		#if typeof(myLoop[i].value)!=1e-139
			try eval(Meta.parse(myLoop[i].name*"="*string(myLoop[i].value)))
			catch
			end
		#end
	end
	prob = ODEProblem(eqDiffToSolve! ,u₀,tspan,p)
	condition(u,t,integrator) = t ∈ dosetimes
	affect!(integrator) = integrator.p[1] = -integrator.p[1]
	cb = DiscreteCallback(condition,affect!)
	sol= DifferentialEquations.solve(prob ,Tsit5(),callback=cb,tstops=dosetimes)
    return sol
end

function PlotSol(sol,myVrArray)
	myGraph=plot(sol,linewidth=2,xaxis="t",label=permutedims(myVrArray),layout=(4,1))
	return myGraph
end

# example d'utilisation du module
#myStrgToSolve1,myLoop,myVrArray=SolveSymbolic(myInputCircuit)
#sol=NumericSolve(myStrgToSolve1,myLoop)
#PlotSol(sol,myVrArray)

end


text:
using AlgebraicMatrices
using Random
using Test

# Set reproducible random number seed
Random.seed!(0)

# Random rationals
function Base.rand(rng::AbstractRNG,
                   ::Random.SamplerType{Rational{T}}) where {T}
    l = lcm(1:10)
    return Rational{T}(T(rand(rng, (-l):l)) // l)
end

@testset "Algebraic matrices as vector spaces D=$D" for D in 0:3
    T = Rational{BigInt}
    niters = 20
    for iter in 1:niters
        print("\r$iter/$niters...")

        sz = rand(0:10, D)
        x = arand(AlgebraicArray{T,D}, sz)
        y = arand(AlgebraicArray{T,D}, sz)
        z = arand(AlgebraicArray{T,D}, sz)
        n = zero(AlgebraicArray{T,D}, sz)
        a = rand(T)
        b = rand(T)

        @test x isa AbstractArray{T,D}
        @test y isa AbstractArray{T,D}
        @test z isa AbstractArray{T,D}
        @test n isa AbstractArray{T,D}
        @test a isa T
        @test b isa T

        @test size(x) == size(y) == size(z) == size(n)

        @test x == x
        @test x == n || x + x ≠ x

        if D == 1
            @test (x < y) + (x == y) + (x > y) == 1
        end

        @test collect(x) isa Array{T,D}
        @test evaluate(x) isa AbstractArray{T,D}
        @test collect(x) == evaluate(x)

        x′ = map(x -> 2x + 1, x)
        @test x′ == map(x -> 2x + 1, evaluate(x))

        @test (x + y) + z == x + (y + z)
        @test n + x == x
        @test x + n == x
        @test x + y == y + x
        @test one(T) * x == x
        @test x * one(T) == x
        @test a * n == n
        @test n * a == n
        if a ≠ zero(T)
            @test inv(a) * (a * x) == x
        end
        @test +x == x
        @test -x == n - x
        @test x + (-x) == n
        @test (-x) + x == n
        @test a * (x + y) == a * x + a * y
        @test (a + b) * x == a * x + b * x

        @test (x == y) == (evaluate(x) == evaluate(y))
        if D == 1
            @test (x < y) == (evaluate(x) < evaluate(y))
            @test (x > y) == (evaluate(x) > evaluate(y))
        end
        @test iszero(evaluate(n))
        @test evaluate(x + y) == evaluate(x) + evaluate(y)
        @test evaluate(a * x) == a * evaluate(x)
        @test evaluate(-x) == -evaluate(x)

        @test simplify(x) == x
    end
    println("\r", " "^40, "\r")
end

@testset "Algebraic matrices as categories D=$D1" for D1 in 1:3
    T = Rational{BigInt}
    niters = 10
    for iter in 1:niters
        print("\r$iter/$niters...")

        D2 = rand(1:3)
        D3 = rand(1:3)
        sz1 = rand(0:10, D1)
        sz2 = (sz1[end], rand(0:10, D2 - 1)...)
        sz3 = (sz2[end], rand(0:10, D3 - 1)...)
        x = arand(AlgebraicArray{T,D1}, sz1)
        x′ = arand(AlgebraicArray{T,D1}, sz1)
        y = arand(AlgebraicArray{T,D2}, sz2)
        z = arand(AlgebraicArray{T,D3}, sz3)
        e1 = one(AlgebraicArray{T,2}, sz1[1], sz1[1])
        e2 = one(AlgebraicArray{T,2}, sz1[end], sz1[end])
        nx = zero(AlgebraicArray{T,D1}, sz1)
        ny = zero(AlgebraicArray{T,D2}, sz2)
        D12 = D1 + D2 - 2
        sz12 = (sz1[1:(end - 1)]..., sz2[2:end]...)
        nxy = zero(AlgebraicArray{T,D12}, sz12)
        a = rand(T)

        if ndims(y) ≥ 2
            @test (x * y) * z == x * (y * z)
        end
        @test e1 * x == x
        @test x * e2 == x
        @test nx * y == nxy
        @test x * ny == nxy
        @test (a * x) * y == a * (x * y)
        @test x * (a * y) == a * (x * y)
        @test (x + x′) * y == x * y + x′ * y

        @test isone(evaluate(e1))

        x′ = evaluate(x)
        y′ = evaluate(y)
        sxy = (size(x)[1:(end - 1)]..., size(y)[2:end]...)
        xy′ = zeros(T, sxy)
        for i in CartesianIndices(size(x)[1:(end - 1)])
            for j in CartesianIndices(size(y)[2:end])
                for k in 1:size(y, 1)
                    xy′[i, j] += x′[i, k] * y′[k, j]
                end
            end
        end
        # @test evaluate(x * y) == evaluate(x) * evaluate(y)
        @test evaluate(x * y) == xy′
    end
    println("\r", " "^40, "\r")
end


text:
using ReplMaker
using .Plots

function _ngspice_parser(str)
    if str == "init" init()
    elseif occursin("plot", str)
        _plotswitch(str)
    else
        cmd(str)
    end
end

function _repl_parser(str)
    if occursin(".sp", str) && str[1:7] == "source "
        source_sp(str[8:end])
    else
        _ngspice_parser(str)
    end
end

interactive() = (init();
           initrepl(_repl_parser,
                  prompt_text  = "NgSpice> ",
                  prompt_color = :yellow,
                  start_key    = '~',
                  mode_name    = "for NgSpice is",
                  ))

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
using TimeSeries, MarketData, Base.Dates
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("find methods") do

    context("find returns correct row numbers array") do
        @fact cl[find(cl .> op)].timestamp[1] --> Date(2000,1,3)
        @fact length(find(cl .> op))          --> 244
    end

    context("findwhen returns correct Dates array") do
       @fact findwhen(cl .> op)[2]      --> Date(2000,1,5)
       @fact length(findwhen(cl .> op)) --> 244
    end
end

facts("split date operations") do

    context("from and to correctly subset non-zero and zero-length time arrays") do
        @fact length(from(cl, Date(2001,12,28))) --> 2
        @fact length(from(cl, Date(2002,1,1)))   --> 0
        @fact length(from(from(cl, Date(2002,1,1)), Date(2012,1,1))) --> 0

        @fact length(to(cl, Date(2000,1,4)))     --> 2
        @fact length(to(cl, Date(1999,1,4)))     --> 0
        @fact length(to(to(cl, Date(1999,1,4)), Date(1912,1,1))) --> 0
    end

    context("when method correctly subset") do
        @fact when(cl, day, 4).timestamp[1]              --> Date(2000,1,4)
        @fact when(cl, dayname, "Friday").timestamp[1]   --> Date(2000,1,7)
        @fact when(cl, week, 5).timestamp[1]             --> Date(2000,1,31)
        @fact when(cl, month, 5).timestamp[1]            --> Date(2000,5,1)
        @fact when(cl, monthname, "June").timestamp[1]   --> Date(2000,6,1)
        @fact when(cl, year, 2001).timestamp[1]          --> Date(2001,1,2)
        @fact when(cl, dayofweek, 1).timestamp[1]        --> Date(2000,1,3)
        # all the days in the nth week of each month
        @fact when(cl, dayofweekofmonth, 5).timestamp[1] --> Date(2000,1,31)
        @fact when(cl, dayofyear, 365).timestamp[1]      --> Date(2001,12,31)
        @fact when(cl, quarterofyear, 4).timestamp[1]    --> Date(2000,10,2)
        @fact when(cl, dayofquarter, 1).timestamp[1]     --> Date(2001,10,1)
    end
end

facts("element wrappers") do

    context("type element wrappers isolate elements") do
        @fact isa(timestamp(cl), Array{Date,1})       --> true
        @fact isa(values(cl), Array{Float64,1})       --> true
        @fact isa(values(ohlc), Array{Float64,2})     --> true
        @fact isa(colnames(cl), Array{String, 1}) --> true
    end
end

facts("head and tail methods") do

    context("head and tail methods work with default n value on single column TimeArray") do
        @fact length(head(cl))      --> 1
        @fact length(head(cl,1))    --> 1
        @fact head(cl).timestamp[1] --> Date(2000,1,3)
        @fact head(cl).values[1]    --> 111.94
        @fact head(cl).meta         --> "AAPL"

        @fact length(tail(cl))      --> 1
        @fact length(tail(cl,1))    --> 1
        @fact tail(cl).timestamp[1] --> Date(2001,12,31)
        @fact tail(cl).values[1]    --> 21.9
        @fact tail(cl).meta         --> "AAPL"
    end

    context("head and tail methods work with default n value on multi column TimeArray") do
        @fact length(head(ohlc))      --> 1
        @fact length(head(ohlc,1))    --> 1
        @fact head(ohlc).timestamp[1] --> Date(2000,1,3)
        @fact head(ohlc).values       --> [104.88  112.5  101.69  111.94]
        @fact head(ohlc).meta         --> "AAPL"

        @fact length(tail(ohlc))      --> 1
        @fact length(tail(ohlc,1))    --> 1
        @fact tail(ohlc).timestamp[1] --> Date(2001,12,31)
        @fact tail(ohlc).values       --> [22.51  22.66  21.83  21.9]
        @fact tail(ohlc).meta         --> "AAPL"
    end

    context("head and tail methods work with custom periods on single column TimeArray") do
        @fact length(head(cl, 2))   --> 2
        @fact length(head(cl, 500)) --> length(cl)
        @fact length(tail(cl, 2))   --> 2
        @fact length(tail(cl, 500)) --> length(cl)
    end

    context("head and tail methods work with custom periods on multi column TimeArray") do
        @fact length(head(ohlc, 2))   --> 2
        @fact length(head(ohlc, 500)) --> length(ohlc)
        @fact length(tail(ohlc, 2))   --> 2
        @fact length(tail(ohlc, 500)) --> length(ohlc)
    end
end


text:
"""
    me_solve_time_independent(rho_in::Array{T1,2},
                                   H::Array{T2,2},
                                   Gamma::Array{T3,3},
                                   rates::Array{T4,1},
                                   t0::AbstractFloat,
                                   tf::AbstractFloat;
                                   tstep::AbstractFloat=0.,
                                   tols::Vector{Float64}=[1e-6,1e-3],
                                   alg = Tsit5()) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: AbstractFloat}

Time independent master equation solver using a non-vectorized algorithm.

## args
REQUIRED
* rho_in:   d x d array, the density matrix of the initial state
* H:        d x d array, system Hamiltonian
* Gamma:    d x d x K array, all K Lindblad operators
* rates:    K x 1 array, dissipative rate for each Lindblad operator
* t0:       float, start time
* tf:       float, end time
KEYWORD OPTIONAL
* tstep:    float, time steps at which output data should be saved
* tols:     2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
* alg:      function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
* iter_num: number of iterations for an adpative solver

## returns
* tvec:     vector of time points where the density matrix has been simulated
* rho_out:  vector of simulated density matricies
"""
function me_solve_time_independent(rho_in::Array{T1,2},
                                   H::Array{T2,2},
                                   Gamma::Array{T3,3},
                                   rates::Array{T4,1},
                                   t0::AbstractFloat,
                                   tf::AbstractFloat;
                                   tstep::AbstractFloat=0.,
                                   tols::Vector{Float64}=[1e-6,1e-3],
                                   alg = Tsit5(),
                                   iter_num=1e5) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: AbstractFloat}

    rho_in       = convert(Array{ComplexF64,2},rho_in)
    H            = convert(Array{ComplexF64,2},H)
    Gamma        = convert(Array{ComplexF64,3},Gamma)
    rates        = convert(Array{Float64,1},rates)
    dRho_L       = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))

    scratchA     = similar(rho_in);
    scratchB     = similar(rho_in);

    Gammas       = [Gamma[:,:,ii].*sqrt(rates[ii]) for ii=1:size(Gamma,3)]
    GammaTs      = [collect(G') for G in Gammas]
    GammaSqs     = [0.5.*G'G for G in Gammas]

    dif_f(du,u,p,t) = dRho!(du,u,H,Gammas,GammaTs,GammaSqs,dRho_L,scratchA,scratchB) # In place
    tspan = (t0,tf)

    prob = ODEProblem{true}(dif_f,rho_in,tspan)

    if tstep < 1e-10
        tstep = (tf-t0)./100
    end

    sol = solve(prob, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], maxiters=iter_num)

    tvec = sol.t
    rho_out = sol.u

    return tvec, rho_out
end
"""
    me_solve_H_time_dependent(rho_in::Array{T1,2},
                                   H::Function,
                                   Gamma::Array{T3,3},
                                   rates::Array{T4,1},
                                   t0::AbstractFloat,
                                   tf::AbstractFloat;
                                   tstep::AbstractFloat=0.,
                                   tols::Vector{Float64}=[1e-6,1e-3],
                                   alg = Tsit5(),
                                   iter_num=1e5,
                                   stop_points=[]
                                   adapt::Bool = true,
                                   δt::Float64 = tols[1],
                                   override::Bool = false) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: AbstractFloat}

Time dependent (Hamiltonian) master equation solver using a non-vectorized algorithm.

## args
REQUIRED
* rho_in:       d x d array, the density matrix of the initial state
* H:            function, system Hamiltonian, takes input of time and returns the
                    Hamiltonian at the given time in matrix form
* Gamma:        d x d x K array, all K Lindblad operators
* rates:        K x 1 array, dissipative rate for each Lindblad operator
* t0:           float, start time
* tf:           float, end time
KEYWORD OPTIONAL
* tstep:        float, time steps at which output data should be saved
* tols:         2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
* alg:          function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
* iter_num:     number of iterations for an adpative solver
* stop_points:  vector of time-points the solver must step through (for adaptive solvers)
* adapt:        Boolean. If false, the solver assumes it is using a fixed
                timestep integration method.
* δt:           Timestep for fixed timestep integration method.
* override:     Booelan. If true, overrides the error checking that the data
                save timestep (tstep) is not less than the fixed integration
                timestep (δt).

## returns
* tvec:         vector of time points where the density matrix has been simulated
* rho_out:      vector of simulated density matricies
"""
function me_solve_H_time_dependent(rho_in::Array{T1,2},
                                   H::Function,
                                   Gamma::Array{T3,3},
                                   rates::Array{T4,1},
                                   t0::AbstractFloat,
                                   tf::AbstractFloat;
                                   tstep::AbstractFloat=0.,
                                   tols::Vector{Float64}=[1e-6,1e-3],
                                   alg = Tsit5(),
                                   iter_num=1e5,
                                   stop_points=[],
                                   adapt::Bool = true,
                                   δt::Float64 = tols[1],
                                   override::Bool = false) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: AbstractFloat}

    rho_in = convert(Array{ComplexF64,2},rho_in)
    Gamma = convert(Array{ComplexF64,3},Gamma)
    rates = convert(Array{Float64,1},rates)

    dRho_L = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))
    H_temp = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))

    scratchA     = similar(rho_in);
    scratchB     = similar(rho_in);

    Gammas       = [Gamma[:,:,ii].*sqrt(rates[ii]) for ii=1:size(Gamma,3)]
    GammaTs      = [collect(G') for G in Gammas]
    GammaSqs     = [0.5.*G'G for G in Gammas]

    dif_f(du,u,p,t) = dRho!(du,u,H,H_temp,Gammas,GammaTs,GammaSqs,dRho_L,scratchA,scratchB,t) # In place
    tspan = (t0,tf)

    prob = ODEProblem{true}(dif_f,rho_in,tspan)

    if tstep < 1e-10
        tstep = (tf-t0)./100
    end

    if adapt
        sol = solve(prob, alg, saveat = tstep, tstops = stop_points, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], maxiters=iter_num)
    else
        if tstep < δt && override == false
            error("You are using a fixed-timestep solver and the save timestep is smaller than solver timestep. This is probably a bad idea. If you insist, set override = true in the input arguments to override this error checking.")
        end
        sol = solve(prob, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], adaptive = adapt, dt = δt)
    end

    tvec = sol.t
    rho_out = sol.u

    return tvec, rho_out
end

"""
    me_solve_H_time_dependent(rho_in::Array{T1,2},
                                   Hops::Array{T5,3},
                                   Hfuncs::Function,
                                   Gamma::Array{T3,3},
                                   rates::Array{T4,1},
                                   t0::AbstractFloat,
                                   tf::AbstractFloat;
                                   tstep::AbstractFloat=0.,
                                   tols::Vector{Float64}=[1e-6,1e-3],
                                   alg = Tsit5(),iter_num=1e5,
                                   stop_points=[],
                                   adapt::Bool = true,
                                   δt::Float64 = tols[1],
                                   override::Bool = false) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: AbstractFloat, T5 <: Number}

Time dependent (Hamiltonian) master equation solver using a non-vectorized algorithm. Instead of inputing the Hamiltonian as a matrix function, it is input as a basis of matrices, and a function that describes the time evolution of the scalar prefactor of each basis element.

This solver also allows for fixed timestep integration.

Input arguments (different from above):
## args
REQUIRED
* Hops:     Array of matrices that describe the Hamiltonian. In the most
                general case is a full basis for operator space.
* Hfuncs:   Function that takes as input a Float (the time) and returns
                a Vector of Floats. Each element of the output vector is
                the value of the scalar prefactor for the corresponding
                basis element in the decomposition of the Hamiltonian.

KEYWORD OPTIONAL
* adapt:    Boolean. If false, the solver assumes it is using a fixed
                timestep integration method.
* δt:       Timestep for fixed timestep integration method.
* override: Booelan. If true, overrides the error checking that the data
                save timestep (tstep) is not less than the fixed integration
                timestep (δt).

## returns
* tvec:         vector of time points where the density matrix has been simulated
* rho_out:      vector of simulated density matricies
"""
function me_solve_H_time_dependent(rho_in::Array{T1,2},
                                   Hops::Array{T5,3},
                                   Hfuncs::Function,
                                   Gamma::Array{T3,3},
                                   rates::Array{T4,1},
                                   t0::AbstractFloat,
                                   tf::AbstractFloat;
                                   tstep::AbstractFloat=0.,
                                   tols::Vector{Float64}=[1e-6,1e-3],
                                   alg = Tsit5(),
                                   iter_num=1e5,
                                   stop_points=[],
                                   adapt::Bool = true,
                                   δt::Float64 = tols[1],
                                   override::Bool = false) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: AbstractFloat, T5 <: Number}

    rho_in = convert(Array{ComplexF64,2},rho_in)
    Gamma = convert(Array{ComplexF64,3},Gamma)
    rates = convert(Array{Float64,1},rates)
    Hops = convert(Array{ComplexF64,3},Hops)
    dRho_L = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))
    H_temp = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))
    Hf_temp = zeros(ComplexF64,size(Hops,3))

    scratchA     = similar(rho_in);
    scratchB     = similar(rho_in);

    Gammas       = [Gamma[:,:,ii].*sqrt(rates[ii]) for ii=1:size(Gamma,3)]
    GammaTs      = [collect(G') for G in Gammas]
    GammaSqs     = [0.5.*G'G for G in Gammas]

    dif_f(du,u,p,t) = dRho!(du,u,Hops,Hfuncs,H_temp,Hf_temp,Gammas,GammaTs,GammaSqs,dRho_L,scratchA,scratchB,t) # In place
    tspan = (t0,tf)

    prob = ODEProblem{true}(dif_f,rho_in,tspan)

    if tstep < 1e-10
        tstep = (tf-t0)./100
    end

    if adapt
        sol = solve(prob, alg, saveat = tstep, tstops = stop_points, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], maxiters=iter_num)
    else
        if tstep < δt && override == false
            error("You are using a fixed-timestep solver and the save timestep is smaller than solver timestep. This is probably a bad idea. If you insist, set override = true in the input arguments to override this error checking.")
        end
        sol = solve(prob, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], adaptive = adapt, dt = δt)
    end

    tvec = sol.t
    rho_out = sol.u

    return tvec, rho_out
end

"""
    me_solve_L_time_dependent(rho_in::Array{T1,2},
                                   H::Array{T2,2},
                                   Gamma::Array{T3,3},
                                   rates::Function,
                                   t0::AbstractFloat,
                                   tf::AbstractFloat;
                                   tstep::AbstractFloat=0.,
                                   tols::Vector{Float64}=[1e-6,1e-3],
                                   alg = Tsit5(),
                                   iter_num=1e5,
                                   stop_points=[],
                                   adapt::Bool = true,
                                   δt::Float64 = tols[1],
                                   override::Bool = false) where {T1 <: Number, T2 <: Number, T3 <: Number}

Time dependent (dissipative rates) master equation solver using a non-vectorized algorithm.

## args
* rho_in:       d x d array, the density matrix of the initial state
* H:            d x d array, system Hamiltonian
* Gamma:        d x d x K array, all K Lindblad operators
* rates:        K x 1 array, dissipative rate for each Lindblad operator, takes input time and returns a vector of the dissipative rates at the given time
* t0:           float, start time
* tf:           float, end time
KEYWORD OPTIONAL
* tstep:        float, time steps at which output data should be saved
* tols:         2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
* alg:          function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
* iter_num:     number of iterations for an adpative solver
* stop_points:  vector of time-points the solver must step through (for adaptive solvers)
* adapt:        Boolean. If false, the solver assumes it is using a fixed
                timestep integration method.
* δt:           Timestep for fixed timestep integration method.
* override:     Booelan. If true, overrides the error checking that the data
                save timestep (tstep) is not less than the fixed integration
                timestep (δt).

## returns
* tvec:         vector of time points where the density matrix has been simulated
* rho_out:      vector of simulated density matricies
"""
function me_solve_L_time_dependent(rho_in::Array{T1,2},
                                   H::Array{T2,2},
                                   Gamma::Array{T3,3},
                                   rates::Function,
                                   t0::AbstractFloat,
                                   tf::AbstractFloat;
                                   tstep::AbstractFloat=0.,
                                   tols::Vector{Float64}=[1e-6,1e-3],
                                   alg = Tsit5(),
                                   iter_num=1e5,
                                   stop_points=[],
                                   adapt::Bool = true,
                                   δt::Float64 = tols[1],
                                   override::Bool = false) where {T1 <: Number, T2 <: Number, T3 <: Number}

    rho_in = convert(Array{ComplexF64,2},rho_in)
    H = convert(Array{ComplexF64,2},H)
    Gamma = convert(Array{ComplexF64,3},Gamma)

    scratchA     = similar(rho_in);
    scratchB     = similar(rho_in);

    Gammas       = [Gamma[:,:,ii] for ii=1:size(Gamma,3)]
    GammaTs      = [collect(G') for G in Gammas]
    GammaSqs     = [0.5.*G'G for G in Gammas]

    dRho_L = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))
    rates_temp = Array{Float64}(undef,1)

    dif_f(du,u,p,t) = dRho!(du,u,H,rates,rates_temp,Gammas,GammaTs,GammaSqs,dRho_L,scratchA,scratchB,t) # In place
    tspan = (t0,tf)

    prob = ODEProblem{true}(dif_f,rho_in,tspan)

    if tstep < 1e-10
        tstep = (tf-t0)./100
    end

    if adapt
        sol = solve(prob, alg, saveat = tstep, tstops = stop_points, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], maxiters=iter_num)
    else
        if tstep < δt && override == false
            error("You are using a fixed-timestep solver and the save timestep is smaller than solver timestep. This is probably a bad idea. If you insist, set override = true in the input arguments to override this error checking.")
        end
        sol = solve(prob, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], adaptive = adapt, dt = δt)
    end

    tvec = sol.t
    rho_out = sol.u

    return tvec, rho_out
end

"""
    me_solve_full_time_dependent(rho_in::Array{T1,2},
                                      H::Function,
                                      Gamma::Array{T3,3},
                                      rates::Function,
                                      t0::AbstractFloat,
                                      tf::AbstractFloat;
                                      tstep::AbstractFloat=0.,
                                      tols::Vector{Float64}=[1e-6,1e-3],
                                      alg = Tsit5(),
                                      stop_points=[],
                                      adapt::Bool = true,
                                      δt::Float64 = tols[1],
                                      override::Bool = false) where {T1 <: Number, T2 <: Number, T3 <: Number}
Time dependent (Hamiltonian and dissipative rates) master equation solver using a non-vectorized algorithm.

## args
REQUIRED
* rho_in:       d x d array, the density matrix of the initial state
* H:            function, system Hamiltonian, takes input of time and returns
                    the Hamiltonian at the given time in matrix form
* Gamma:        d x d x K array, all K Lindblad operators
* rates:        K x 1 array, dissipative rate for each Lindblad operator,
                    takes input time and returns a vector of the dissipative
                    rates at the given time
* t0:           float, start time
* tf:           float, end time
KEYWORD OPTIONAL
* tstep:        float, time steps at which output data should be saved
* tols:         2 x 1 array, vector of solver tolernaces in the
                    order [abstol, reltol]
* alg:          function, algorithm from DifferentialEquations for the solver
                    to use, default is Tsit5
* iter_num:     number of iterations for an adpative solver
* stop_points:  vector of time-points the solver must step through
                    (for adaptive solvers)
* adapt:        Boolean. If false, the solver assumes it is using a fixed
                timestep integration method.
* δt:           Timestep for fixed timestep integration method.
* override:     Booelan. If true, overrides the error checking that the data
                save timestep (tstep) is not less than the fixed integration
                timestep (δt).

## returns
* tvec:         vector of time points where the density matrix has been simulated
* rho_out:      vector of simulated density matricies
"""
function me_solve_full_time_dependent(rho_in::Array{T1,2},
                                      H::Function,
                                      Gamma::Array{T3,3},
                                      rates::Function,
                                      t0::AbstractFloat,
                                      tf::AbstractFloat;
                                      tstep::AbstractFloat=0.,
                                      tols::Vector{Float64}=[1e-6,1e-3],
                                      alg = Tsit5(),
                                      iter_num=1e5,
                                      stop_points=[],
                                      adapt::Bool = true,
                                      δt::Float64 = tols[1],
                                      override::Bool = false) where {T1 <: Number, T2 <: Number, T3 <: Number}

    rho_in = convert(Array{ComplexF64,2},rho_in)
    Gamma = convert(Array{ComplexF64,3},Gamma)

    dRho_L = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))
    H_temp = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))
    rates_temp = Array{Float64}(undef,1)

    scratchA     = similar(rho_in);
    scratchB     = similar(rho_in);

    Gammas       = [Gamma[:,:,ii] for ii=1:size(Gamma,3)]
    GammaTs      = [collect(G') for G in Gammas]
    GammaSqs     = [0.5.*G'G for G in Gammas]

    dif_f(du,u,p,t) = dRho!(du,u,H,H_temp,rates,rates_temp,Gammas,GammaTs,GammaSqs,dRho_L,scratchA,scratchB,t) # In place
    tspan = (t0,tf)

    prob = ODEProblem{true}(dif_f,rho_in,tspan)

    if tstep < 1e-10
        tstep = (tf-t0)./100
    end

    if adapt
        sol = solve(prob, alg, saveat = tstep, tstops = stop_points, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], maxiters=iter_num)
    else
        if tstep < δt && override == false
            error("You are using a fixed-timestep solver and the save timestep is smaller than solver timestep. This is probably a bad idea. If you insist, set override = true in the input arguments to override this error checking.")
        end
        sol = solve(prob, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2], adaptive = adapt, dt = δt)
    end

    tvec = sol.t
    rho_out = sol.u

    return tvec, rho_out
end

"""
    me_solve_time_independent_vec(rho_in::Array{T1,2},
                                       H::Array{T2,2},
                                       Gamma::Array{T3,3},
                                       rates::Array{T4,1},
                                       t0::AbstractFloat,
                                       tf::AbstractFloat;
                                       tstep::AbstractFloat=0.,
                                       tols::Vector{Float64}=[1e-6,1e-3],
                                       alg = Tsit5()) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: AbstractFloat}

Time independent master equation solver using a non-vectorized algorithm but with vector ODE solvers

## args
REQUIRED
* rho_in:       d x d array, the density matrix of the initial state
* H:            d x d array, system Hamiltonian
* Gamma:        d x d x K array, all K Lindblad operators
* rates:        K x 1 array, dissipative rate for each Lindblad operator
* t0:           float, start time
* tf:           float, end time
KEYWORD OPTIONAL
* tstep:        float, time steps at which output data should be saved
* tols:         2 x 1 array, vector of solver tolernaces in the
                    order [abstol, reltol]
* alg:          function, algorithm from DifferentialEquations for the
                    solver to use, default is Tsit5

## returns
* tvec:         vector of time points where the density matrix has been simulated
* rho_out:      vector of simulated density matricies
"""
# LUKE: I think we can remove this functionality.
# function me_solve_time_independent_vec(rho_in::Array{T1,2},
#                                        H::Array{T2,2},
#                                        Gamma::Array{T3,3},
#                                        rates::Array{T4,1},
#                                        t0::AbstractFloat,
#                                        tf::AbstractFloat;
#                                        tstep::AbstractFloat=0.,
#                                        tols::Vector{Float64}=[1e-6,1e-3],
#                                        alg = Tsit5()) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: AbstractFloat}
#
#     rho_in = convert(Array{ComplexF64,2},rho_in)
#     H = convert(Array{ComplexF64,2},H)
#     Gamma = convert(Array{ComplexF64,3},Gamma)
#     rates = convert(Array{Float64,1},rates)
#
#     rho_in_vec = rho_in[:]
#
#     dRho_L = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))
#     rho_temp = Array{ComplexF64}(undef,size(rho_in,1),size(rho_in,2))
#
#     dif_f(du,u,p,t) = dRho_vec(du,u,H,Gamma,rates,dRho_L,rho_temp) # In place
#     tspan = (t0,tf)
#
#     prob = ODEProblem{true}(dif_f,rho_in_vec,tspan)
#
#     if tstep < 1e-10
#         tstep = (tf-t0)./100
#     end
#
#     sol = solve(prob, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2])
#
#     tvec = sol.t
#     rho_out_vec = sol.u
#
#     rho_out = Array{Array{ComplexF64,2}}(undef,length(tvec))
#     for jj = 1:1:length(tvec)
#         rho_out[jj] = reshape(rho_out_vec[jj][:],size(H))
#     end
#
#     return tvec, rho_out
# end

# """
# Time dependent solver for Gaussian continuous variable systems using the covariance matrix approach. Defined in terms of the canonical X and P quadratures, the covariance matrix elements are defined by
#   C_ab = <O_aO_b> - <O_a><O_b> where vec(O) = [X,P] with
#   X = (a + a^†)/√2, P = -i(a - a^†)/√2 with a the lowering operator of a harmonic oscilaltor.
# Input Parameters:
# REQUIRED
#   C_in: Array 2 x 2, initial condition for the covariance matrix
#   a_in: vector 2 x 1, initial condition for the average value vector [<X>,<P>]
#   h: Function that returnas as output a matrix 2 x 2, Hamiltonian evolution term for the covariance matrix. For H = ω a^†a + (λa^2 + conj(λ)a^†^2)/2, h = [Im(λ), (ω - Re(λ);-(ω + Re(λ)), -Im(λ)].
#       Can use h = create_h_C(ω(t),λ(t),t) to generate the correct h.
#   drv: Function that returns as output a vector 2 x 1, linear drive term. For H = αa + conj(α)a^†, drv = √2[Im(α), Re(α)]
#   t0: float, start time
#   tf: float, end time
# KEYWORD OPTIONAL
#   tstep: float, time steps at which output data should be saved
#   tols: 2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
#   alg: function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
# """
#
# function CV_solve_full_time_dependent(C_in::Array{T,2},a_in::Vector{Float64},h::Function,drv::Function,t0::AbstractFloat,tf::AbstractFloat; tstep::AbstractFloat=0.,tols::Vector{Float64}=[1e-6,1e-3],alg = Tsit5()) where {T <: Number}
#
#   C_in = convert(Array{ComplexF64,2}, C_in)
#   if tstep < 1e-10
#         tstep = (tf-t0)./100
#     end
#   tspan = (t0,tf)
#
#   # Covariance matrix
#   h_tempL = Array{Float64}(undef,size(C_in,1),size(C_in,2))
#   # h_tempR = Array{Float64}(undef,size(C_in,1),size(C_in,2))
#     dif_C(du,u,p,t) = dCV(du,u,h,t,h_tempL) # In place
#
#     prob_C = ODEProblem{true}(dif_C,C_in,tspan)
#
#     sol_C = solve(prob_C, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2])
#
#   # Average
#   h_temp = Array{Float64}(undef,size(C_in,1),size(C_in,2))
#   drv_temp = Array{Float64}(undef,size(a_in,1))
#   dif_a(du,u,p,t) = d_av(du,u,h,t,drv,h_temp,drv_temp) # In place
#
#     prob_a = ODEProblem{true}(dif_a,a_in,tspan)
#
#     sol_a = solve(prob_a, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2])
#
# return sol_C.t, sol_C.u, sol_a.t, sol_a.u
# end


text:
# This function creates bipolar (i.e. signed) log of the input.

function bilog(datain, dbrange=40)
    datain = real.(datain)
    maxabs = abs.(maximum(datain))
    minabs = abs.(minimum(datain))

    if minabs > maxabs
        maxabs = minabs
    end

    norm_mag = (abs.(datain) .+ eps()) / maxabs
    out = sign.(datain) .* (clamp.(20 * log10.(norm_mag), -dbrange, 0) .+ dbrange)

    return out
end


text:
using DataFrames
using HDF5

function main()
    table_path = ARGS[1]
    id = h5read(table_path,"id")
    x = h5read(table_path,"x")
    y = h5read(table_path,"y")

    df = DataFrame()
    df[:id] = id
    df[:x] = x
    df[:y] = y
    A = Array(Float64,size(x,1))
    tic()
    for i in 2:size(x,1)-1
        A[i] = (df[:x][i-1]+2*df[:x][i]+df[:x][i+1])/4.
    end
    t = toq()
    println("DataFrames.jl WMA took: ", t, " val", sum(A))
end

main()

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
cd(@__DIR__); include("setups/path128.jl");
gr(dpi = 200)

## (a)
Uf = unitary_folding_operator(W, GP; ϵ = 0.38, J = 2)' # ϵa ≈ 0.38⋅∥𝛟₁∥_{∞}
plt = plot(size = (500, 400))
    heatmap!(Uf, c = :viridis, ratio = 1, yflip = true, xlim = [1, N], clims = (-1, 1))
    xticks!([1; 16:16:128], vcat(string(1), [string(k) for k in 16:16:128]))
    yticks!([1; 16:16:128], vcat(string(1), [string(k) for k in 16:16:128]))
savefig(plt, "../figs/Path128_Uf_j2.png")

## (b)
Uf = unitary_folding_operator(W, GP; ϵ = 0.38, J = 3)'
plt = plot(size = (500, 400))
    heatmap!(Uf, c = :viridis, ratio = 1, yflip = true, xlim = [1, N], clims = (-1, 1))
    xticks!([1; 16:16:128], vcat(string(1), [string(k) for k in 16:16:128]))
    yticks!([1; 16:16:128], vcat(string(1), [string(k) for k in 16:16:128]))
savefig(plt, "../figs/Path128_Uf_j3.png")


text:
module CoreMessage

using ..Ahorn, Maple

# Render both vanilla and everest Core Messages, but only allow placement of Everest version as it allows custom dialog
const placements = Ahorn.PlacementDict(
    "Core Message (Everest)" => Ahorn.EntityPlacement(
        Maple.EverestCoreMessage
    )
)

coreMessageUnion = Union{Maple.EverestCoreMessage, Maple.CoreMessage}

function Ahorn.selection(entity::coreMessageUnion)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::coreMessageUnion, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.speechBubble, -12, -12)

end


text:
# This file contains types and structures for technical DEA model
"""
    AbstractTechnicaDEAlModel
An abstract type representing a technical DEA model.
"""
abstract type AbstractTechnicalDEAModel  <: AbstractDEAModel end

"""
    efficiency(model::AbstractTechnicalDEAModel)
Return efficiency scores of a technical DEA model.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> efficiency(deaio)
11-element Vector{Float64}:
 1.0
 0.6222896790980051
 0.8198562443845464
 1.0
 ⋮
 0.7576690895651103
 0.8201058201058201
 0.49056603773584917
 1.0
```
"""
efficiency(model::AbstractTechnicalDEAModel) = model.eff

"""
    slacks(model::AbstractTechnicalDEAModel, slack::Symbol)
Return slacks of a technical DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> slacks(deaio, :X)
11×2 Matrix{Float64}:
  0.0          0.0
 -4.41868e-15  0.0
  0.0          8.17926e-15
 -8.03397e-16  0.0
  ⋮
  1.60679e-15  0.0
  1.64021      0.0
  9.68683e-15  0.0
  0.0          4.0

julia> slacks(deaio, :Y)
11×1 Matrix{Float64}:
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
```
"""
function slacks(model::AbstractTechnicalDEAModel, slack::Symbol)::Matrix

    if slack == :X
        return model.slackX
    elseif slack == :Y
        return model.slackY
    end

    throw(ArgumentError("`slack` must be :X or :Y"));

end

"""
    targets(model::AbstractTechnicalDEAModel, target::Symbol)
Return targets of a technical DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> targets(deaio, :X)
11×2 Matrix{Float64}:
  5.0      13.0
  9.95663   7.46748
 13.1177   21.3163
 17.0      15.0
  ⋮
 20.4571   16.6687
 28.7037   11.4815
 20.6038   12.2642
  5.0      13.0

julia> targets(deaio, :Y)
11×1 Matrix{Float64}:
 12.0
 14.0
 25.0
 26.0
  ⋮
 30.0
 31.0
 26.0
 12.0
```
"""
function targets(model::AbstractTechnicalDEAModel, target::Symbol)::Matrix

    if target == :X
        return model.Xtarget
    elseif target == :Y
        return model.Ytarget
    end

    throw(ArgumentError("`target` must be :X or :Y"));

end


text:
"""
    power_spectrum_3d(::AtmosLESConfigType, u, v, w, L, dim, nor)

Calculates the Powerspectrum from the 3D velocity fields `u`, `v`, `w`. Inputs
need to be equi-spaced and the domain is assumed to be the same size and
have the same number of points in all directions.

# Arguments
- L: size domain
- dim: number of points
- nor: normalization factor
"""
function power_spectrum_3d(::AtmosLESConfigType, u, v, w, L, dim, nor)
    u_fft = fft(u)
    v_fft = fft(v)
    w_fft = fft(w)

    mu = Array((abs.(u_fft) / size(u, 1)^3))
    mv = Array((abs.(v_fft) / size(v, 1)^3))
    mw = Array((abs.(w_fft) / size(w, 1)^3))
    if mod(dim, 2) == 0
        rx = range(0, stop = dim - 1, step = 1) .- dim / 2 .+ 1
        ry = range(0, stop = dim - 1, step = 1) .- dim / 2 .+ 1
        rz = range(0, stop = dim - 1, step = 1) .- dim / 2 .+ 1
        R_x = circshift(rx', (1, dim / 2 + 1))
        R_y = circshift(ry', (1, dim / 2 + 1))
        R_z = circshift(rz', (1, dim / 2 + 1))
        k_nyq = dim / 2
    else
        rx = range(0, stop = dim - 1, step = 1) .- (dim - 1) / 2
        ry = range(0, stop = dim - 1, step = 1) .- (dim - 1) / 2
        rz = range(0, stop = dim - 1, step = 1) .- (dim - 1) / 2
        R_x = circshift(rx', (1, (dim + 1) / 2))
        R_y = circshift(ry', (1, (dim + 1) / 2))
        R_z = circshift(rz', (1, (dim + 1) / 2))
        k_nyq = (dim - 1) / 2
    end
    r = zeros(size(rx, 1), size(ry, 1), size(rz, 1))
    for i in 1:size(rx, 1), j in 1:size(ry, 1), l in 1:size(rz, 1)
        r[i, j, l] = sqrt(R_x[i]^2 + R_y[j]^2 + R_z[l]^2)
    end
    dx = 2 * pi / L
    k = range(1, stop = k_nyq, step = 1) .* dx
    endk = size(k, 1)
    contribution = zeros(size(k, 1))
    spectrum = zeros(size(k, 1))
    for N in 2:Int64(k_nyq - 1)
        for i in 1:size(rx, 1), j in 1:size(ry, 1), l in 1:size(rz, 1)
            if (r[i, j, l] * dx <= (k'[N + 1] + k'[N]) / 2) &&
               (r[i, j, l] * dx > (k'[N] + k'[N - 1]) / 2)
                spectrum[N] =
                    spectrum[N] + mu[i, j, l]^2 + mv[i, j, l]^2 + mw[i, j, l]^2
                contribution[N] = contribution[N] + 1
            end
        end
    end
    for i in 1:size(rx, 1), j in 1:size(ry, 1), l in 1:size(rz, 1)
        if (r[i, j, l] * dx <= (k'[2] + k'[1]) / 2)
            spectrum[1] =
                spectrum[1] + mu[i, j, l]^2 + mv[i, j, l]^2 + mw[i, j, l]^2
            contribution[1] = contribution[1] + 1
        end
    end
    for i in 1:size(rx, 1), j in 1:size(ry, 1), l in 1:size(rz, 1)
        if (r[i, j, l] * dx <= k'[endk]) &&
           (r[i, j, l] * dx > (k'[endk] + k'[endk - 1]) / 2)
            spectrum[endk] =
                spectrum[endk] + mu[i, j, l]^2 + mv[i, j, l]^2 + mw[i, j, l]^2
            contribution[endk] = contribution[endk] + 1
        end
    end
    spectrum = spectrum .* 2 .* pi .* k .^ 2 ./ (contribution .* dx .^ 3) ./ nor

    return spectrum, k
end