text:
import LanguageServer.Document
import LanguageServer.get_text
import LanguageServer.update
import LanguageServer.get_offset
import LanguageServer.get_line_offsets
import LanguageServer.get_position_at
import LanguageServer.get_open_in_editor
import LanguageServer.set_open_in_editor
import LanguageServer.is_workspace_file

s1 = """
123456
abcde
ABCDEFG
"""
d1 = Document("untitled", s1, false)
@test get_text(d1) == s1
@test get_offset(d1, 1, 4) == 4
@test get_offset(d1, 2, 2) == 9
@test get_line_offsets(d1) == [1, 8, 14, 21 + 1]
@test get_position_at(d1, 1) == (1, 1)
@test get_position_at(d1, 8) == (2, 1)
@test get_position_at(d1, 15) == (3, 2)
@test get_open_in_editor(d1) == false
set_open_in_editor(d1, true)
@test get_open_in_editor(d1) == true
set_open_in_editor(d1, false)
@test get_open_in_editor(d1) == false
@test is_workspace_file(d1) == false

s2 = """
12μ456
abηde
ABCDEFG"""
d2 = Document("untitled", s2, true)
@test get_offset(d2, 1, 4) == 5
@test get_offset(d2, 2, 2) == 10
@test get_line_offsets(d2) == [1, 9, 16]
@test get_line_offsets(d2) == [nextind(d2._content,0,1), nextind(d2._content,0,8), nextind(d2._content,0,14)]
@test get_position_at(d2, 1) == (1, 1)
@test get_position_at(d2, 9) == (2, 1)
@test get_position_at(d2, 17) == (3, 2)
@test is_workspace_file(d2) == true

update(d2, 2, 2, 0, "12")
@test get_line_offsets(d2) == [nextind(d2._content, 0, 1), nextind(d2._content, 0, 8), nextind(d2._content, 0, 16)]

update(d2, 1, 3, 3, "abcdef")
@test get_line_offsets(d2) == [nextind(d2._content, 0, 1), nextind(d2._content, 0, 11), nextind(d2._content, 0, 19)]

update(d2, 2, 5, 4, "xyz")
@test get_line_offsets(d2) == [nextind(d2._content, 0, 1), nextind(d2._content, 0, 11)]

update(d2, 1, 1, 0, "PRE")
@test get_text(d2) == "PRE12abcdef6\na12bxyzABCDEFG"

update(d2, 2, 15, 0, "POST")
@test get_text(d2) == "PRE12abcdef6\na12bxyzABCDEFGPOST"

s4 = "1234\r\nabcd"
d4 = Document("untitled", s4, false)
@test get_line_offsets(d4) == [1, 7]

s5 = "1234\nabcd\n"
d5 = Document("untitled", s5, false)
@test get_line_offsets(d5) == [1, 6, 11]

s6 = "\n"
d6 = Document("untitled", s6, false)
@test get_line_offsets(d6) == [1, 2]


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

JLLWrappers.@generate_main_file_header("SCALAPACK32")
JLLWrappers.@generate_main_file("SCALAPACK32", UUID("aabda75e-bfe4-5a37-92e3-ffe54af3c273"))
end  # module SCALAPACK32_jll

ext: jl


text:
GetH(x) = (suff(x) == BigFloat) ? convert(BigFloat,10^(-precision(BigFloat)/10)) : 1e-6

# """
#     suff(x) -> Type
# If `x` stores BigFloats, `suff` returns BigFloat, else `suff` returns `Float64`.
# """
# suff(x::BigFloat) = BigFloat
# suff(x::Float32) = Float32
# suff(x::Float16) = Float16
# suff(x::Real) = Float64
# suff(x::Num) = Num
# suff(x::Complex) = suff(real(x))
# suff(x::AbstractArray) = suff(x[1])
# suff(x::DataFrame) = suff(x[1,1])
# suff(x::Tuple) = suff(x...)
# suff(args...) = try suff(promote(args...)[1]) catch;  suff(args[1]) end
# # Allow for differentiation through suff arrays.
# suff(x::ForwardDiff.Dual) = typeof(x)
# suff(x::ReverseDiff.TrackedReal) = typeof(x)

floatify(x::AbstractArray{<:AbstractFloat}) = x;   floatify(x::AbstractArray) = float.(x)
floatify(x::AbstractFloat) = x;                     floatify(x::Number) = float(x)
floatify(x) = float.(x)

"""
    Unpack(Z::AbstractVector{S}) where S <: Union{AbstractVector,Tuple} -> Matrix
Converts vector of vectors to a matrix whose n-th column corresponds to the n-th component of the inner vectors.
"""
@inline function Unpack(Z::AbstractVector{S}) where S <: Union{AbstractVector{<:Number},Tuple}
    N = length(Z);      M = length(Z[1])
    A = Array{suff(Z)}(undef,N,M)
    @inbounds for i in Base.OneTo(N)
        for j in Base.OneTo(M)
            A[i,j] = Z[i][j]
        end
    end;    A
end
Unpack(Z::AbstractVector{<:Number}) = Z

Unwind(M::AbstractMatrix{<:Number}) = Unwind(collect(eachrow(M)))
Unwind(X::AbstractVector{<:AbstractVector{<:Number}}) = reduce(vcat, X)
Unwind(X::AbstractVector{<:Number}) = X

Windup(X::AbstractVector{<:Number}, n::Int) = n < 2 ? X : [X[(1+(i-1)*n):(i*n)] for i in 1:Int(length(X)/n)]

ToCols(M::Matrix) = Tuple(M[:,i] for i in 1:size(M,2))

ValToBool(x::Val{true}) = true
ValToBool(x::Val{false}) = false

function GetMethod(tol::Real)
    if tol > 1e-8
        Tsit5()
    elseif tol < 1e-11
        Vern9()
    else
        Vern7()
    end
end

# Check for length
PromoteStatic(X::AbstractArray, inplace::Bool=true) = length(X) > 90 ? X : PromoteStatic(X, Val(inplace))

# No checking for length
PromoteStatic(X::AbstractArray, mutable::Val{true}) = _PromoteMutable(X)
PromoteStatic(X::AbstractArray, mutable::Val{false}) = _PromoteStatic(X)

_PromoteMutable(X::AbstractVector, Length=length(X)) = MVector{Length}(X)
_PromoteMutable(X::AbstractArray, Size=size(X)) = MArray{Tuple{Size...}}(X)
_PromoteStatic(X::AbstractVector, Length=length(X)) = SVector{Length}(X)
_PromoteStatic(X::AbstractArray, Size=size(X)) = SArray{Tuple{Size...}}(X)

# Surely, this can be made more efficient?
SplitAfter(n::Int) = X->(X[1:n], X[n+1:end])

"""
    invert(F::Function, x::Number; tol::Real=GetH(x)) -> Real
Finds ``z`` such that ``F(z) = x`` to a tolerance of `tol` for continuous ``F`` using Roots.jl. Ideally, `F` should be monotone and there should only be one correct result.
"""
function invert(F::Function, x::Number, Domain::Tuple{<:Number,<:Number}=(zero(suff(x)), 1e4*one(suff(x)));
                    tol::Real=GetH(x), meth::Roots.AbstractUnivariateZeroMethod=Roots.Order1())
    @assert Domain[1] < Domain[2]
    try
        if meth isa Roots.AbstractNonBracketing
            find_zero(z-> F(z) - x, 0.5one(suff(x)), meth; xatol=tol)
        else
            find_zero(z-> F(z) - x, Domain, meth; xatol=tol)
        end
    catch err
        @warn "invert() errored: $(nameof(typeof(err))). Assuming result is bracketed by $Domain and falling back to Bisection-like method."
        find_zero(z-> F(z) - x, Domain, Roots.AlefeldPotraShi(); xatol=tol)
    end
end
# function invert(F::Function, x::Number; tol::Real=GetH(x)*100, meth::Roots.AbstractUnivariateZeroMethod=Order1())
#     find_zero(z-> F(z) - x, 0.8one(suff(x)), meth; xatol=tol)
# end

"""
    ConfAlpha(n::Real)
Probability volume outside of a confidence interval of level n⋅σ where σ is the standard deviation of a normal distribution.
"""
ConfAlpha(n::Real) = 1 - ConfVol(n)

"""
    ConfVol(n::Real)
Probability volume contained in a confidence interval of level n⋅σ where σ is the standard deviation of a normal distribution.
"""
function ConfVol(n::Real)
    if abs(n) ≤ 8
        return erf(n / sqrt(2))
    else
        println("ConfVol: Float64 precision not enough for n = $n. Returning BigFloat instead.")
        return ConfVol(BigFloat(n))
    end
end
ConfVol(n::BigFloat) = erf(n / sqrt(BigFloat(2)))

InvConfVol(q::Real; kwargs...) = sqrt(2) * erfinv(q)
InvConfVol(x::BigFloat; tol::Real=GetH(x)) = invert(ConfVol, x; tol=tol)

ChisqCDF(k::Int, x::Real) = gamma_inc(k/2., x/2., 0)[1]
# ChisqCDF(k::Int, x::Real) = cdf(Chisq(k), x)
ChisqCDF(k::Int, x::BigFloat) = gamma_inc(BigFloat(k)/2., x/2., 0)[1]

InvChisqCDF(k::Int, p::Real; kwargs...) = 2gamma_inc_inv(k/2., p, 1-p)
InvChisqCDF(k::Int, p::BigFloat; tol::Real=GetH(p)) = invert(x->ChisqCDF(k, x), p; tol=tol)

InnerProduct(Mat::AbstractMatrix, Y::AbstractVector) = transpose(Y) * Mat * Y
# InnerProduct(Mat::PDMats.PDMat, Y::AbstractVector) = (R = Mat.chol.U * Y;  dot(R,R))

import Base.==
==(DS1::DataSet, DS2::DataSet) = xdata(DS1) == xdata(DS2) && ydata(DS1) == ydata(DS2) && ysigma(DS1) == ysigma(DS2)
==(DS1::DataSetExact, DS2::DataSet) = DS2 == DS1
function ==(DS1::DataSet, DS2::DataSetExact)
    if !(xdist(DS2) isa InformationGeometry.Dirac && ydist(DS2) isa MvNormal)
        return false
    elseif xdata(DS1) == xdata(DS2) && ydata(DS1) == ydata(DS2) && yInvCov(DS1) == yInvCov(DS2)
        return true
    else
        false
    end
end
==(DS1::AbstractDataSet, DS2::AbstractDataSet) = xdist(DS1) == xdist(DS2) && ydist(DS1) == ydist(DS2)

"""
    GetArgSize(model::ModelOrFunction; max::Int=100)
Returns tuple `(xdim,pdim)` associated with the method `model(x,p)`.
"""
GetArgSize(model::Function; max::Int=100) = GetArgSize(model, Val(MaximalNumberOfArguments(model) > 2); max=max)
function GetArgSize(model::Function, inplace::Val{false}; max::Int=100)
    try return (1, GetArgLength(p->model(rand(),p); max=max)) catch; end
    for i in 2:max
        try return (i,GetArgLength(p->model(rand(i),p); max=max)) catch; end
    end;    throw("Wasn't able to find config for max=$max.")
end
function GetArgSize(model!::Function, inplace::Val{true}; max::Int=100)
    try return (1, GetArgLength((Res,p)->model!(Res,rand(),p); max=max)) catch; end
    for i in 2:max
        try return (i,GetArgLength((Res,p)->model!(Res,rand(i),p); max=max)) catch; end
    end;    throw("Wasn't able to find config for max=$max.")
end
GetArgSize(model::ModelMap; max::Int=100) = (model.xyp[1], model.xyp[3])

normalize(x::AbstractVector{<:Number}, scaling::Float64=1.0) = (scaling / norm(x)) * x
function normalizeVF(u::AbstractVector{<:Number}, v::AbstractVector{<:Number}, scaling::Float64=1.0)
    newu = u;    newv = v
    for i in 1:length(u)
        factor = sqrt(u[i]^2 + v[i]^2)
        newu[i] = (scaling/factor)*u[i]
        newv[i] = (scaling/factor)*v[i]
    end
    newu, newv
end
function normalizeVF(u::AbstractVector{<:Number},v::AbstractVector{<:Number},PlanarCube::HyperCube,scaling::Float64=1.0)
    length(PlanarCube) != 2 && throw("normalizeVF: Cube not planar.")
    newu = u;    newv = v
    Widths = CubeWidths(PlanarCube) |> normalize
    for i in 1:length(u)
        factor = sqrt(u[i]^2 + v[i]^2)
        newu[i] = (scaling/factor)*u[i] * Widths[1]
        newv[i] = (scaling/factor)*v[i] * Widths[2]
    end
    newu, newv
end

"""
    BlockMatrix(M::AbstractMatrix, N::Int)
Returns matrix which contains `N` many blocks of the matrix `M` along its diagonal.
"""
function BlockMatrix(M::AbstractMatrix, N::Int)
    Res = zeros(size(M,1)*N,size(M,2)*N)
    for i in 1:N
        Res[((i-1)*size(M,1) + 1):(i*size(M,1)),((i-1)*size(M,1) + 1):(i*size(M,1))] = M
    end;    Res
end
BlockMatrix(M::Diagonal, N::Int) = Diagonal(repeat(M.diag, N))

"""
    BlockMatrix(A::AbstractMatrix, B::AbstractMatrix)
Constructs blockdiagonal matrix from `A` and `B`.
"""
function BlockMatrix(A::AbstractMatrix, B::AbstractMatrix)
    Res = zeros(suff(A), size(A,1)+size(B,1), size(A,2)+size(B,2))
    Res[1:size(A,1), 1:size(A,1)] = A
    Res[size(A,1)+1:end, size(A,1)+1:end] = B
    Res
end
BlockMatrix(A::Diagonal, B::Diagonal) = Diagonal(vcat(A.diag, B.diag))

BlockMatrix(A::AbstractMatrix, B::AbstractMatrix, args...) = BlockMatrix(BlockMatrix(A,B), args...)


text:
module TestMeasures

using Fairness, Test
using MLJBase, MLJModels

include("calcmetrics.jl")
include("boolmetrics.jl")
include("fairness.jl")
include("metricWrapper.jl")

end  # module


text:
export Diag,
       DiagTensor

# Diag can have either Vector storage, in which case
# it is a general Diag tensor, or scalar storage,
# in which case the diagonal has a uniform value
struct Diag{ElT,VecT} <: TensorStorage{ElT}
  data::VecT
  function Diag{ElT,VecT}(data) where {ElT,VecT<:AbstractVector{ElT}}
    return new{ElT,VecT}(data)
  end
  function Diag{ElT,ElT}(data) where {ElT}
    return new{ElT,ElT}(data)
  end
end

Diag(data::VecT) where {VecT<:AbstractVector{ElT}} where {ElT} = Diag{ElT,VecT}(data)

Diag(data::ElT) where {ElT<:Number} = Diag{ElT,ElT}(data)

function Diag{ElR}(data::AbstractVector{ElT}) where {ElR<:Number,ElT<:Number}
  ElT == ElR ? Diag(data) : Diag(ElR.(data))
end

Diag(::Type{ElT},
     n::Integer) where {ElT<:Number} = Diag(zeros(ElT,n))

Diag(x::ElT,
     n::Integer) where {ElT<:Number} = Diag(fill(x,n))

Base.copy(D::Diag) = Diag(copy(data(D)))

const NonuniformDiag{ElT,VecT} = Diag{ElT,VecT} where {VecT<:AbstractVector}

const UniformDiag{ElT,VecT} = Diag{ElT,VecT} where {VecT<:Number}

Base.getindex(D::UniformDiag,i::Int) = data(D)

Base.setindex!(D::UniformDiag,val,i::Int) = error("Cannot set elements of a uniform Diag storage")

Base.complex(::Type{Diag{ElT,Vector{ElT}}}) where {ElT} = Diag{complex(ElT),Vector{complex(ElT)}}
Base.complex(::Type{Diag{ElT,ElT}}) where {ElT} = Diag{complex(ElT),complex(ElT)}

# Deal with uniform Diag conversion
Base.convert(::Type{<:Diag{ElT,VecT}},D::Diag) where {ElT,VecT} = Diag(convert(VecT,data(D)))

# TODO: write in terms of ::Int, not inds
Base.similar(D::NonuniformDiag) = Diag(similar(data(D)))
#Base.similar(D::NonuniformDiag,inds) = Diag(similar(data(D),minimum(dims(inds))))
#function Base.similar(D::Type{<:NonuniformDiag{ElT,VecT}},inds) where {ElT,VecT}
#  return Diag(similar(VecT,diaglength(inds)))
#end

Base.similar(D::UniformDiag{ElT}) where {ElT} = Diag(zero(ElT))
Base.similar(D::UniformDiag,inds) = similar(D)
Base.similar(::Type{<:UniformDiag{ElT}},inds) where {ElT} = Diag(zero(ElT))

Base.similar(D::Diag,n::Int) = Diag(similar(data(D),n))

Base.similar(D::Diag,::Type{ElR},n::Int) where {ElR} = Diag(similar(data(D),ElR,n))

# TODO: make this work for other storage besides Vector
Base.zeros(::Type{<:NonuniformDiag{ElT}},dim::Int64) where {ElT} = Diag(zeros(ElT,dim))
Base.zeros(::Type{<:UniformDiag{ElT}},dim::Int64) where {ElT} = Diag(zero(ElT))

Base.:*(D::Diag,x::Number) = Diag(x*data(D))
Base.:*(x::Number,D::Diag) = D*x

#
# Type promotions involving Diag
# Useful for knowing how conversions should work when adding and contracting
#

function Base.promote_rule(::Type{<:UniformDiag{ElT1}},
                           ::Type{<:UniformDiag{ElT2}}) where {ElT1,ElT2}
  ElR = promote_type(ElT1,ElT2)
  return Diag{ElR,ElR}
end

function Base.promote_rule(::Type{<:NonuniformDiag{ElT1,VecT1}},
                           ::Type{<:NonuniformDiag{ElT2,VecT2}}) where {ElT1,VecT1<:AbstractVector,
                                                                        ElT2,VecT2<:AbstractVector}
  ElR = promote_type(ElT1,ElT2)
  VecR = promote_type(VecT1,VecT2)
  return Diag{ElR,VecR}
end

# This is an internal definition, is there a more general way?
#Base.promote_type(::Type{Vector{ElT1}},
#                  ::Type{ElT2}) where {ElT1<:Number,
#                                       ElT2<:Number} = Vector{promote_type(ElT1,ElT2)}
#
#Base.promote_type(::Type{ElT1},
#                  ::Type{Vector{ElT2}}) where {ElT1<:Number,
#                                               ElT2<:Number} = promote_type(Vector{ElT2},ElT1)

# TODO: how do we make this work more generally for T2<:AbstractVector{S2}?
# Make a similar_type(AbstractVector{S2},T1) -> AbstractVector{T1} function?
function Base.promote_rule(::Type{<:UniformDiag{ElT1,VecT1}},
                           ::Type{<:NonuniformDiag{ElT2,Vector{ElT2}}}) where {ElT1,VecT1<:Number,
                                                                               ElT2}
  ElR = promote_type(ElT1,ElT2)
  VecR = Vector{ElR}
  return Diag{ElR,VecR}
end

function Base.promote_rule(::Type{DenseT1},
                           ::Type{<:NonuniformDiag{ElT2,VecT2}}) where {DenseT1<:Dense,
                                                                        ElT2,VecT2<:AbstractVector}
  return promote_type(DenseT1,Dense{ElT2,VecT2})
end

function Base.promote_rule(::Type{DenseT1},
                           ::Type{<:UniformDiag{ElT2,VecT2}}) where {DenseT1<:Dense,
                                                                     ElT2,VecT2<:Number}
  return promote_type(DenseT1,ElT2)
end

# Convert a Diag storage type to the closest Dense storage type
dense(::Type{<:NonuniformDiag{ElT,VecT}}) where {ElT,VecT} = Dense{ElT,VecT}
dense(::Type{<:UniformDiag{ElT}}) where {ElT} = Dense{ElT,Vector{ElT}}

const DiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:Diag}
const NonuniformDiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where
                                               {StoreT<:NonuniformDiag}
const UniformDiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where
                                               {StoreT<:UniformDiag}

Base.IndexStyle(::Type{<:DiagTensor}) = IndexCartesian()

# TODO: this needs to be better (promote element type, check order compatibility,
# etc.
function Base.convert(::Type{<:DenseTensor{ElT,N}}, T::DiagTensor{ElT,N}) where {ElT<:Number,N}
  return dense(T)
end

# These are rules for determining the output of a pairwise contraction of NDTensors
# (given the indices of the output tensors)
function contraction_output_type(TensorT1::Type{<:DiagTensor},
                                 TensorT2::Type{<:DenseTensor},
                                 IndsR::Type)
  return similar_type(promote_type(TensorT1,TensorT2),IndsR)
end
contraction_output_type(TensorT1::Type{<:DenseTensor},
                        TensorT2::Type{<:DiagTensor},
                        IndsR::Type) = contraction_output_type(TensorT2,TensorT1,IndsR)

# This performs the logic that DiagTensor*DiagTensor -> DiagTensor if it is not an outer
# product but -> DenseTensor if it is
# TODO: if the tensors are both order 2 (or less), or if there is an Index replacement,
# then they remain diagonal. Should we limit DiagTensor*DiagTensor to cases that
# result in a DiagTensor, for efficiency and type stability? What about a general
# SparseTensor result?
function contraction_output_type(TensorT1::Type{<:DiagTensor{<:Number,N1}},
                                 TensorT2::Type{<:DiagTensor{<:Number,N2}},
                                 IndsR::Type) where {N1,N2}
  if ValLength(IndsR)===Val{N1+N2}
    # Turn into is_outer(inds1,inds2,indsR) function?
    # How does type inference work with arithmatic of compile time values?
    return similar_type(dense(promote_type(TensorT1,TensorT2)),IndsR)
  end
  return similar_type(promote_type(TensorT1,TensorT2),IndsR)
end

# The output must be initialized as zero since it is sparse, cannot be undefined
contraction_output(T1::DiagTensor,T2::Tensor,indsR) = zero_contraction_output(T1,T2,indsR)
contraction_output(T1::Tensor,T2::DiagTensor,indsR) = contraction_output(T2,T1,indsR)

function contraction_output(T1::DiagTensor,
                            T2::DiagTensor,
                            indsR)
  return zero_contraction_output(T1,T2,indsR)
end

function array(T::DiagTensor{ElT,N}) where {ElT,N}
  return array(dense(T))
end
matrix(T::DiagTensor{<:Number,2}) = array(T)
vector(T::DiagTensor{<:Number,1}) = array(T)

function Base.Array{ElT,N}(T::DiagTensor{ElT,N}) where {ElT,N}
  return array(T)
end

function Base.Array(T::DiagTensor{ElT,N}) where {ElT,N}
  return Array{ElT,N}(T)
end

# Needed to get slice of DiagTensor like T[1:3,1:3]
function Base.similar(T::DiagTensor{<:Number,N},
                      ::Type{ElR},
                      inds::Dims{N}) where {ElR<:Number,N}
  return tensor(similar(store(T),ElR,minimum(inds)),inds)
end

"""
getdiagindex(T::DiagTensor,i::Int)

Get the ith value along the diagonal of the tensor.
"""
getdiagindex(T::DiagTensor{<:Number},ind::Int) = store(T)[ind]

"""
setdiagindex!(T::DiagTensor,i::Int)

Set the ith value along the diagonal of the tensor.
"""
setdiagindex!(T::DiagTensor{<:Number},val,ind::Int) = (store(T)[ind] = val)

"""
setdiag(T::UniformDiagTensor,val)

Set the entire diagonal of a uniform DiagTensor.
"""
setdiag(T::UniformDiagTensor,val) = tensor(Diag(val),inds(T))

Base.@propagate_inbounds function Base.getindex(T::DiagTensor{ElT,N},
                                                inds::Vararg{Int,N}) where {ElT,N}
  if all(==(inds[1]),inds)
    return getdiagindex(T,inds[1])
  else
    return zero(eltype(ElT))
  end
end
Base.@propagate_inbounds Base.getindex(T::DiagTensor{<:Number,1},ind::Int) = store(T)[ind]
Base.@propagate_inbounds Base.getindex(T::DiagTensor{<:Number,0}) = store(T)[1]

# Set diagonal elements
# Throw error for off-diagonal
Base.@propagate_inbounds function Base.setindex!(T::DiagTensor{<:Number,N},
                                                 val,inds::Vararg{Int,N}) where {N}
  all(==(inds[1]),inds) || error("Cannot set off-diagonal element of Diag storage")
  setdiagindex!(T,val,inds[1])
  return T
end
Base.@propagate_inbounds Base.setindex!(T::DiagTensor{<:Number,1},val,ind::Int) = ( store(T)[ind] = val )
Base.@propagate_inbounds Base.setindex!(T::DiagTensor{<:Number,0},val) = ( store(T)[1] = val )

function Base.setindex!(T::UniformDiagTensor{<:Number,N},val,inds::Vararg{Int,N}) where {N}
  error("Cannot set elements of a uniform Diag storage")
end

# TODO: make a fill!! that works for uniform and non-uniform
#Base.fill!(T::DiagTensor,v) = fill!(store(T),v)

function dense(::Type{<:Tensor{ElT,N,StoreT,IndsT}}) where {ElT,N,
                                                            StoreT<:Diag,IndsT}
  return Tensor{ElT,N,dense(StoreT),IndsT}
end

# convert to Dense
function dense(T::TensorT) where {TensorT<:DiagTensor}
  R = zeros(dense(TensorT),inds(T))
  for i = 1:diaglength(T)
    setdiagindex!(R,getdiagindex(T,i),i)
  end
  return R
end

function outer!(R::DenseTensor{<:Number,NR},
                T1::DiagTensor{<:Number,N1},
                T2::DiagTensor{<:Number,N2}) where {NR,N1,N2}
  for i1 = 1:diaglength(T1), i2 = 1:diaglength(T2)
    indsR = CartesianIndex{NR}(ntuple(r -> r ≤ N1 ? i1 : i2, Val(NR)))
    R[indsR] = getdiagindex(T1,i1)*getdiagindex(T2,i2)
  end
  return R
end

# TODO: write an optimized version of this?
function outer!(R::DenseTensor{ElR},
                T1::DenseTensor,
                T2::DiagTensor) where {ElR}
  R .= zero(ElR)
  outer!(R,T1,dense(T2))
  return R
end

function outer!(R::DenseTensor{ElR},
                T1::DiagTensor,
                T2::DenseTensor) where {ElR}
  R .= zero(ElR)
  outer!(R,dense(T1),T2)
  return R
end

# Right an in-place version
function outer(T1::DiagTensor{ElT1,N1},
               T2::DiagTensor{ElT2,N2}) where {ElT1,ElT2,N1,N2}
  indsR = unioninds(inds(T1),inds(T2))
  R = tensor(Dense(zeros(promote_type(ElT1,ElT2),dim(indsR))),indsR)
  outer!(R,T1,T2)
  return R
end

function Base.permutedims!(R::DiagTensor{<:Number,N},
                           T::DiagTensor{<:Number,N},
                           perm::NTuple{N,Int},f::Function=(r,t)->t) where {N}
  # TODO: check that inds(R)==permute(inds(T),perm)?
  for i=1:diaglength(R)
    @inbounds setdiagindex!(R,f(getdiagindex(R,i),getdiagindex(T,i)),i)
  end
  return R
end

function Base.permutedims(T::DiagTensor{<:Number,N},
                          perm::NTuple{N,Int},f::Function=identity) where {N}
  R = similar(T,permute(inds(T),perm))
  permutedims!(R,T,perm,f)
  return R
end

function Base.permutedims(T::UniformDiagTensor{ElT,N},
                          perm::NTuple{N,Int},
                          f::Function=identity) where {ElR,ElT,N}
  R = tensor(Diag(f(getdiagindex(T,1))),permute(inds(T),perm))
  return R
end

# Version that may overwrite in-place or may return the result
function permutedims!!(R::NonuniformDiagTensor{<:Number,N},
                       T::NonuniformDiagTensor{<:Number,N},
                       perm::NTuple{N,Int},
                       f::Function=(r,t)->t) where {N}
  permutedims!(R,T,perm,f)
  return R
end

function permutedims!!(R::UniformDiagTensor{ElR,N},
                       T::UniformDiagTensor{ElT,N},
                       perm::NTuple{N,Int},
                       f::Function=(r,t)->t) where {ElR,ElT,N}
  R = tensor(Diag(f(getdiagindex(R,1),getdiagindex(T,1))),inds(R))
  return R
end

function Base.permutedims!(R::DenseTensor{ElR,N},
                           T::DiagTensor{ElT,N},
                           perm::NTuple{N,Int},
                           f::Function = (r,t)->t) where {ElR,ElT,N}
  for i = 1:diaglength(T)
    @inbounds setdiagindex!(R,f(getdiagindex(R,i),getdiagindex(T,i)),i)
  end
  return R
end

function permutedims!!(R::DenseTensor{ElR,N},
                       T::DiagTensor{ElT,N},
                       perm::NTuple{N,Int},f::Function=(r,t)->t) where {ElR,ElT,N}
  permutedims!(R,T,perm,f)
  return R
end

function _contract!!(R::UniformDiagTensor{ElR,NR},labelsR,
                     T1::UniformDiagTensor{<:Number,N1},labelsT1,
                     T2::UniformDiagTensor{<:Number,N2},labelsT2) where {ElR,NR,N1,N2}
  if NR==0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    R = setdiag(R,diaglength(T1)*getdiagindex(T1,1)*getdiagindex(T2,1))
  else
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    R = setdiag(R,getdiagindex(T1,1)*getdiagindex(T2,1))
  end
  return R
end

function contract!(R::DiagTensor{ElR,NR},labelsR,
                   T1::DiagTensor{<:Number,N1},labelsT1,
                   T2::DiagTensor{<:Number,N2},labelsT2) where {ElR,NR,N1,N2}
  if NR==0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    Rdiag = zero(ElR)
    for i = 1:diaglength(T1)
      Rdiag += getdiagindex(T1,i)*getdiagindex(T2,i)
    end
    setdiagindex!(R,Rdiag,1)
  else
    min_dim = min(diaglength(T1),diaglength(T2))
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    for i = 1:min_dim
      setdiagindex!(R,getdiagindex(T1,i)*getdiagindex(T2,i),i)
    end
  end
  return R
end

function contract!(C::DenseTensor{ElC,NC},Clabels,
                   A::DiagTensor{ElA,NA},Alabels,
                   B::DenseTensor{ElB,NB},Blabels) where {ElA,NA,
                                                          ElB,NB,
                                                          ElC,NC}
  if all(i -> i < 0, Blabels)
    # If all of B is contracted
    # TODO: can also check NC+NB==NA
    min_dim = minimum(dims(B))
    if length(Clabels) == 0
      # all indices are summed over, just add the product of the diagonal
      # elements of A and B
      for i = 1:min_dim
        setdiagindex!(C,getdiagindex(C,1)+getdiagindex(A,i)*getdiagindex(B,i),1)
      end
    else
      # not all indices are summed over, set the diagonals of the result
      # to the product of the diagonals of A and B
      # TODO: should we make this return a Diag storage?
      for i = 1:min_dim
        setdiagindex!(C,getdiagindex(A,i)*getdiagindex(B,i),i)
      end
    end
  else
    astarts = zeros(Int,length(Alabels))
    bstart = 0
    cstart = 0
    b_cstride = 0
    nbu = 0
    for ib = 1:length(Blabels)
      ia = findfirst(==(Blabels[ib]),Alabels)
      if !isnothing(ia)
        b_cstride += Base.stride(B,ib)
        bstart += astarts[ia]*Base.stride(B,ib)
      else
        nbu += 1
      end
    end

    c_cstride = 0
    for ic = 1:length(Clabels)
      ia = findfirst(==(Clabels[ic]),Alabels)
      if !isnothing(ia)
        c_cstride += Base.stride(C,ic)
        cstart += astarts[ia]*Base.stride(C,ic)
      end
    end

    # strides of the uncontracted dimensions of
    # B
    bustride = zeros(Int,nbu)
    custride = zeros(Int,nbu)
    # size of the uncontracted dimensions of
    # B, to be used in CartesianIndices
    busize = zeros(Int,nbu)
    n = 1
    for ib = 1:length(Blabels)
      if Blabels[ib] > 0
        bustride[n] = Base.stride(B,ib)
        busize[n] = size(B,ib)
        ic = findfirst(==(Blabels[ib]),Clabels)
        custride[n] = Base.stride(C,ic)
        n += 1
      end
    end

    boffset_orig = 1-sum(Base.strides(B))
    coffset_orig = 1-sum(Base.strides(C))
    cartesian_inds = CartesianIndices(Tuple(busize))
    for inds in cartesian_inds
      boffset = boffset_orig
      coffset = coffset_orig
      for i in 1:nbu
        ii = inds[i]
        boffset += ii*bustride[i]
        coffset += ii*custride[i]
      end
      for j in 1:diaglength(A)
        C[cstart+j*c_cstride+coffset] += getdiagindex(A,j)*
                                         B[bstart+j*b_cstride+boffset]
      end
    end
  end
end
contract!(C::DenseTensor,Clabels,
          A::DenseTensor,Alabels,
          B::DiagTensor,Blabels) = contract!(C,Clabels,
                                             B,Blabels,
                                             A,Alabels)


text:
include("utils_dg.jl")

K = 6  # Number of elements

n = 4
nn = 1
r = jacobiGL(0, 0, n)
V = vandermonde(r, 0, 0, n)
fr = collect(range(-1,1,length=100))
fV = vandermonde(fr, 0, 0, n)
refine = fV * inv(V)
# cell average
avg = Diagonal(zeros(n+1))
avg[1:nn] .= 1
p0 = V * avg * inv(V)
# polyone
nn = 2
tmp = zeros(n+1)
tmp[1:nn] .= 1
avg = Diagonal(tmp)
p1 = V * avg * inv(V)
# poly two
nn = 2+1
tmp = zeros(n+1)
tmp[1:nn] .= 1
avg = Diagonal(tmp)
p2 = V * avg * inv(V)

# poly 6
Ω = (; a=0, b=2π, periodic = true)
a = Ω.a
b = Ω.b
mesh = Mesh(K, n, Ω.a, Ω.b, periodic = Ω.periodic)
x = mesh.x
weights = sum(mesh.M, dims = 1)
Δx = x[2] - x[1]

# inexact integration
DM = Diagonal(sum(mesh.M, dims = 1)[:])
mesh.M .= DM
mesh.Mi .= inv(DM)
mesh.lift[:,1] .= mesh.Mi[1,:]
mesh.lift[:,end] .= mesh.Mi[end,:]

b = Ω.b
a = Ω.a

u_first  = mesh.D *  @.sin( 2 * 2π/(b-a) * x)
using Random
Random.seed!(123)
u_first = compute_volume_terms(u_first, mesh) + 0.1 * compute_surface_terms(randn(size(u_first)), mesh)
uo = refine * u_first
x = refine * x

DVu = refine * compute_volume_terms(u_first, mesh)
DAu = refine * compute_surface_terms(u_first, mesh)
Du = DVu + DAu

fig = Figure(resolution = (1700+600, 1000+400))

lims = (extrema(x)..., extrema(uo)...)
axo = Axis(fig[1,1], xlabel = L"x", ylabel = L"y", ylabelsize = 22,
    xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
    xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10)

axo.limits = (extrema(x)..., (extrema(uo) .* 1.1)...)
axo.title = "Original"
axo.titlesize = 32
for i in 1:K
    lines!(axo, x[:,i], uo[:,i], linewidth = 5)
end

ax0 = Axis(fig[1,2], xlabel = L"x", ylabel = L"y", ylabelsize = 22,
    xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
    xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10, ylims = (-1.1, 1.1))
lims2 = (extrema(x)..., (extrema(DVu) .* 1.1)...)
ax0.limits = lims2
ax0.title = "Volume Contribution"
ax0.titlesize = 32
for i in 1:K
    lines!(ax0, x[:,i], DVu[:,i], linewidth = 5)
end

ax1 = Axis(fig[2,1], xlabel = L"x", ylabel = L"y", ylabelsize = 22,
    xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
    xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10, ylims = (-1.1, 1.1))
ax1.limits = (extrema(x)..., (extrema(DAu) .* 1.1)...)
ax1.title = "Surface Contribution"
ax1.titlesize = 32
for i in 1:K
    lines!(ax1, x[:,i], DAu[:,i], linewidth = 5)
end

ax3 = Axis(fig[2,2], xlabel = L"x", ylabel = L"y", ylabelsize = 22,
    xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
    xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10, ylims = (-1.1, 1.1))
ax3.limits = (extrema(x)..., (extrema(Du) .* 1.1)...)
ax3.title = "Derivative"
ax3.titlesize = 32
for i in 1:K
    lines!(ax3, x[:,i], Du[:,i], linewidth = 5)
end

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end
using DataStructures
const IntSet = DataStructures.IntSet
import Compat: String
using Primes

tests = ["int_set",
         "deque",
         "circ_deque",
         "sorted_containers",
         "stack_and_queue",
         "accumulator",
         "classified_collections",
         "disjoint_set",
         "binheap",
         "mutable_binheap",
         "default_dict",
         "ordered_dict",
         "ordered_set",
         "trie",
         "list",
         "multi_dict",
         "circular_buffer",
         "sorting",
         "priorityqueue"]

if length(ARGS) > 0
    tests = ARGS
end

for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    println("$fp ...")
    include(fp)
end


text:
using Test
import Roots.fzero

## Test `fzero` interface to `find_zero`
## test `fzeros` interface for functions

@testset "fzero(s) interface" begin

    ## test interface
    @test fzero(sin, 3)  ≈ pi # order0
    @test fzero(sin, (3,3.1), order=1)  ≈ pi # use order if specified
    @test fzero(sin, (3,4))  ≈ pi  # bracketing
    @test fzero(sin, 3, 4)  ≈ pi  # bracketing
    @test fzero(sin, [3,4])  ≈ pi  # bracketing
    @test_throws ArgumentError fzero(sin, (3,3.1)) # not a bracket
    @test fzero(sin, cos, 3)  ≈ pi # newton

    ## order keyword:
    for o in  keys(Roots._method_lookup)
        @test fzero(x -> x^3 - x, 0.9, order=o)  ≈ 1.0
    end

    ## bypass order keyword
    for M in [Roots.Order0(), Roots.Order1(), Roots.Order1B(), Roots.Order2(), Roots.Order2B()]
        @test fzero(x -> x^3 - x, 0.7, M)  ≈ 1.0
    end

    for M in [Roots.Order1(), Roots.Order1B(), Roots.Order2(), Roots.Order2B()]
        N = Roots.Bisection()
        @test fzero(x -> x^3 - x, 0.7, M, N)  ≈ 1.0
    end

    ### test tolerances
    fn, xstar, x0, br = x -> x^5 - x - 1, 1.1673039782614187, (1.0,1.1), [1.0, 2.0]
    @test fzero(fn, x0, order=1) ≈ xstar
    @test !(fzero(fn, x0, order=1, xatol=1e-4, atol=1) ≈ xstar)
    @test_throws Roots.ConvergenceFailed fzero(fn, x0, order=1, maxevals=3)
    @test !(fzero(fn, x0, order=1, maxevals=3, atol=1e-3) ≈ xstar)
    ## @test !(fzero(fn, x0, order=1, maxfnevals=3, atol=1e-3) ≈ xstar) removed maxfnevals

    ## Infinities
    ## f(b) = Inf
    f =  x -> x + exp(x)
    @test fzero(f, -1e6, 1e6)  ≈ -0.5671432904097838
    ## f(a) = Inf
    f =  x -> 1/x - 1
    @test fzero(f, 0, 2)  ≈ 1.0

    ## test infinite range
    @test fzero(x -> x, -Inf, Inf)  ≈ 0.0

    ##################################################
    ## fzeros function
    ## test interface
    @test fzeros(sin, 3, 4) ≈ [pi]
    @test fzeros(sin, (3, 4)) ≈ [pi]
    @test fzeros(sin, [3, 4]) ≈ [pi]

    @test length(fzeros(x -> exp(x) - x^4, -10, 20)) == 3
    rts = 1:5
    @test all((abs.(fzeros(x -> prod([x-r for r in rts]),0,10)) .- collect(1:5)) .<= 1e-15)

    fn = x -> cos(10*pi*x)
    @test length(fzeros(fn, 0, 1)) == 10

    ### issue with fzeros and roots near 'b'
    @test 0 <  maximum(fzeros(x -> sin(x) - 1/1000*x, 0, pi)) < pi

end


text:
using DelimitedFiles
using Plots, StatsPlots, Measures
gr()

M = 1 # markersize

include("bonds.jl")
# MPO results

numplots = 6

meraplots = Matrix{Any}(undef, (numplots, 3))

data = ["_mera_triv.txt", "_mera_z2.txt", "_mera_u1.txt"]
dims = [meradims_triv, meradims_z2, meradims_u1]
name = ["Trivial", "Z2 symmetry", "U1 symmetry"]
for j = 1:3
    s = data[j]
    data_tensorkit = try
        readdlm("results/tensorkit"*s)
    catch
        nothing
    end
    ymin = minimum(data_tensorkit)
    ymax = maximum(data_tensorkit)
    data_itensors = try
        readdlm("results/itensors"*s) - readdlm("results/itensors"*s[1:end-4]*"_gc.txt")
    catch
        nothing
    end
    if !isnothing(data_itensors)
        ymin = min(ymin, minimum(data_itensors))
        ymax = max(ymax, maximum(data_itensors))
    end
    data_tenpy = try
        readdlm("results/tenpy"*s)
    catch
        nothing
    end
    if !isnothing(data_tenpy)
        ymin = min(ymin, minimum(data_tenpy))
        ymax = max(ymax, maximum(data_tenpy))
    end
    data_tensornetwork = try
        readdlm("results/tensornetwork"*s)
    catch
        nothing
    end
    if !isnothing(data_tensornetwork)
        ymin = min(ymin, minimum(data_tensornetwork))
        ymax = max(ymax, maximum(data_tensornetwork))
    end
    ymin = 10^(floor(log10(ymin)))
    ymax = 10^(ceil(log10(ymax)))

    for i = 1:numplots
        Dmera = dims[j][1][i]

        p = boxplot(data_tensorkit[:,i], label="TensorKit.jl", markersize = M)
        !isnothing(data_itensors) &&
            boxplot!(p, data_itensors[:,i], label="ITensors.jl", markersize = M)
        !isnothing(data_tenpy) &&
            boxplot!(p, data_tenpy[:,i], label="Tenpy", markersize = M)
        !isnothing(data_tensornetwork) &&
            boxplot!(p, data_tensornetwork[:,i], label="TensorNetwork", markersize = M)
        plot!(p, yscale = :log10, yminorticks = true, xticks = [], xgrid = false, xshowaxis = false, ylim = (ymin, ymax))
        if i == 1
            plot!(p, yguide = name[j], yguide_position = :left, yguidefontsize = 10, ytickfontsize = 8, left_margin = 6mm)
            plot!(p, title = "χ = $Dmera", titlefontsize = 10)
        else
            plot!(p, title = "$Dmera", titlefontsize = 10)
            plot!(p, left_margin = -5mm, yformatter = _->"", yshowaxis = false)
        end
        if i == numplots && j == 3
            plot!(p, legendfontsize = 10, legend = :bottomright)
        else
            plot!(p, legend = false)
        end
        meraplots[i, j] = p
    end
end
plotgrid = plot(meraplots..., layout=grid(3,numplots), link = :y, thickness_scaling = 1, size = (1500, 800))
savefig(plotgrid, "meraresults.pdf")


text:
#tpsv
#trsv

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
print("getting started...\r")
using Pkg
Pkg.activate(".")
print("getting started...\r")
using JuDoc
using MLJ
MLJ.color_off();
println(rpad("ready ✓", 50))

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    printv2mol, printv2sdf, sdfilewriter

function printv2atoms(io::IO, mol::SDFile)
    for (i, atom) in enumerate(nodeattrs(mol))
        x, y, z = atom.coords
        xyzsym = @sprintf "%10.4f%10.4f%10.4f %-3s" x y z string(atom.symbol)
        println(io, "$(xyzsym) 0  0  0  0  0  0  0  0  0  0  0  0")
    end
    return
end

function printv2atoms(io::IO, mol::SMILES, coords)
    for (i, atom) in enumerate(nodeattrs(mol))
        x, y = coords[i, 1:2]
        z = 0.0
        xyzsym = @sprintf "%10.4f%10.4f%10.4f %-3s" x y z string(atom.symbol)
        println(io, "$(xyzsym) 0  0  0  0  0  0  0  0  0  0  0  0")
    end
    return
end

function printv2bonds(io::IO, mol::SDFile)
    for (i, (u, v)) in enumerate(edgesiter(mol))
        bond = edgeattr(mol, i)
        uv = @sprintf "%3d%3d%3d%3d  0  0  0" u v bond.order bond.notation
        println(io, uv)
    end
    return
end

function printv2bonds(io::IO, mol::SMILES, styles)
    for (i, (u, v)) in enumerate(edgesiter(mol))
        bond = edgeattr(mol, i)
        if styles[i] in (2, 7)
            f, s = (v, u)
            notation = styles[i] - 1
        else
            f, s = (u, v)
            notation = styles[i]
        end
        uv = @sprintf "%3d%3d%3d%3d  0  0  0" f s bond.order notation
        println(io, uv)
    end
    return
end

function printv2properties(io::IO, mol::GraphMol)
    charges = Tuple{Int,Int}[]
    radicals = Tuple{Int,Int}[]
    masses = Tuple{Int,Float64}[]
    for (i, atom) in enumerate(nodeattrs(mol))
        atom.charge == 0 || push!(charges, (i, atom.charge))
        atom.multiplicity == 1 || push!(radicals, (i, atom.multiplicity))
        atom.mass === nothing || push!(masses, (i, atom.mass))
    end
    if !isempty(charges)
        head = @sprintf "M  CHG%3d" length(charges)
        print(io, head)
        for (i, chg) in charges
            rcd = @sprintf " %3d %3d" i chg
            print(io, rcd)
        end
        println(io)
    end
    if !isempty(radicals)
        head = @sprintf "M  RAD%3d" length(radicals)
        print(io, head)
        for (i, rad) in radicals
            rcd = @sprintf " %3d %3d" i rad
            print(io, rcd)
        end
        println(io)
    end
    if !isempty(masses)
        head = @sprintf "M  ISO%3d" length(masses)
        print(io, head)
        for (i, iso) in masses
            rcd = @sprintf " %3d %3d" i iso
            print(io, rcd)
        end
        println(io)
    end
    return
end

function printv2data(io::IO, mol::GraphMol)
    for (key, val) in mol.attributes
        println(io, "> <$(string(key))>")
        println(io, string(val))
        println(io, "")
    end
    return
end

function printv2mol(io::IO, mol::GraphMol)
    println(io)
    println(io, "MolecularGraph.jl version $(Util.VERSION)")
    println(io)
    ncnt = nodecount(mol)
    ecnt = edgecount(mol)
    header = @sprintf "%3d%3d  0  0  0  0  0  0  0  0999 V2000" ncnt ecnt
    println(io, header)
    if nodeattrtype(mol) === SmilesAtom
        coords, styles = coordgen(mol)
        printv2atoms(io, mol, coords)
        printv2bonds(io, mol, styles)
    else
        printv2atoms(io, mol)
        printv2bonds(io, mol)
    end
    printv2properties(io, mol)
    println(io, "M  END")
    return
end

function printv2mol(mol::GraphMol)
    buf = IOBuffer(write=true)
    printv2mol(buf, mol)
    res = String(take!(buf))
    close(buf)
    return res
end

function printv2sdf(io::IO, mol::GraphMol)
    printv2mol(io, mol)
    printv2data(io, mol)
    println(io, raw"$$$$")
    return
end

"""
    sdfilewriter(io::IO, mols)
    sdfilewriter(filename::AbstractString, mols)

Write molecule data to the output stream as a SDFile format file.
"""
sdfilewriter(io::IO, mols; writer=printv2sdf) = writer.((io,), mols)
sdfilewriter(filename::AbstractString, mols; kwargs...
    ) = sdfilewriter(open(filename, "w"), mols; kwargs...)


text:
struct StackPlay end

mutable struct Stack{T}
    content::Union{T, Nothing}
    link::Union{Id{Stack{T}}, Nothing}
end

Stack{T}() where T = Stack{T}(nothing, nothing)

function hear(s::Scene{Stack{T}}, msg::Tuple{Symbol, Union{Id, T}}) where T
    (type, m) = msg

    if type === :push!
        if !isnothing(my(s).content)
            my(s).link = invite!(s, Stack(my(s).content, nothing))
        end

        my(s).content = m
    elseif type === :pop! && !isnothing(my(s).content)
        say(s, m, my(s).content)
        my(s).content = nothing

        isnothing(my(s).link) || forward!(s, my(s).link)
    else
        error("Can't handle $type")
    end
end

function hear(s::Scene{StackPlay}, ::Genesis!)
    stack = invite!(s, Stack{Symbol}())

    say(s, stack, (:push!, :a))
    @test ask(s, stack, (:pop!, me(s)), Symbol) == :a

    say(s, stack, (:push!, :b))
    @test ask(s, stack, (:pop!, me(s)), Symbol) == :b

    stack = invite!(s, Stack{Int}())

    for i in 1:5
        say(s, stack, (:push!, i))
    end

    @test ask(s, stack, (:pop!, me(s)), Int) == 5

    say(s, stack, (:push!, 6))
    @test ask(s, stack, (:pop!, me(s)), Int) == 6
    @test ask(s, stack, (:pop!, me(s)), Int) == 4

    leave!(s)
end

@testset LuvvyTestSet "Actors Stack" begin
    play!(StackPlay())
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

import Compat.ASCIIString

using GrammarExpts, ACASXProblem
using ExprSearch, DerivationTrees
using RLESUtils, CPUTimeUtils, Loggers, Observers
using DataFrameSets

"""
reps=inner loop reps before taking a time
trials=number of times to take
so the number of evals is trials*rep
"""
function test1(;
    data::AbstractString="dasc",
    n_trials::Int64=50,
    n_reps::Int64=2500,
    verbose::Bool=true,
    b_get_fitness::Bool=false,
    b_emulated_fitness::Bool=false,
    b_predicts::Bool=false,
    b_emulated_predicts::Bool=false)

    problem = ACASXClustering(data)
    grammar = get_grammar(problem)
    p = DerivTreeParams(grammar, 20)
    tree = DerivationTree(p)

    #setup logs
    logs = TaggedDFLogger()
    add_folder!(logs, "params", [ASCIIString, Any], [:name, :value])
    add_folder!(logs, "rand", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "get_expr", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "get_fitness", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "to_function", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "to_function2", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "emulated_fitness", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "predicts", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "emulated_predicts", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "emulated_predicts2", [Int64, Float64], [:trial, :t_elapsed])
    add_folder!(logs, "metrics", [Int64, Float64], [:trial, :t_elapsed])
    observer = Observer()
    add_observer(observer, "params", push!_f(logs, "params"))
    add_observer(observer, "rand", push!_f(logs, "rand"))
    add_observer(observer, "get_expr", push!_f(logs, "get_expr"))
    add_observer(observer, "get_fitness", push!_f(logs, "get_fitness"))
    add_observer(observer, "to_function", push!_f(logs, "to_function"))
    add_observer(observer, "to_function2", push!_f(logs, "to_function2"))
    add_observer(observer, "emulated_fitness", push!_f(logs, "emulated_fitness"))
    add_observer(observer, "predicts", push!_f(logs, "predicts"))
    add_observer(observer, "emulated_predicts", push!_f(logs, "emulated_predicts"))
    add_observer(observer, "emulated_predicts2", push!_f(logs, "emulated_predicts2"))
    add_observer(observer, "metrics", push!_f(logs, "metrics"))

    @notify_observer(observer, "params", ["data", data])
    @notify_observer(observer, "params", ["n_trials", n_trials])
    @notify_observer(observer, "params", ["n_reps", n_reps])

    #time rand!
    println("timing rand!")
    for trial = 1:n_trials
        tstart = CPUtime_start()
        for i = 1:n_reps
            #########################
            rand!(tree, retries=typemax(Int64))
            #########################
        end
        t_elapsed = CPUtime_elapsed_s(tstart)
        @notify_observer(observer, "rand", [trial, t_elapsed])
    end
    #time get_expr
    println("timing get_expr")
    for trial = 1:n_trials
        tstart = CPUtime_start()
        for i = 1:n_reps
            #########################
            expr = get_expr(tree)
            #########################
        end
        t_elapsed = CPUtime_elapsed_s(tstart)
        @notify_observer(observer, "get_expr", [trial, t_elapsed])
    end
    expr = get_expr(tree)

    if b_get_fitness
        #time get_fitness
        println("timing get_fitness")
        for trial = 1:n_trials
            tstart = CPUtime_start()
            for i = 1:n_reps
                #########################
                fitness = get_fitness(problem, expr)
                #########################
            end
            t_elapsed = CPUtime_elapsed_s(tstart)
            @notify_observer(observer, "get_fitness", [trial, t_elapsed])
        end
    end
    #time to_function
    println("timing to_function")
    for trial = 1:n_trials
        tstart = CPUtime_start()
        for i = 1:n_reps
            #########################
            f = to_function(problem, expr)
            #########################
        end
        t_elapsed = CPUtime_elapsed_s(tstart)
        @notify_observer(observer, "to_function", [trial, t_elapsed])
    end
    #time to_function2
    D_1 = getrecords(problem.Dl, 1) #only the first record
    println("timing to_function2")
    for trial = 1:n_trials
        tstart = CPUtime_start()
        for i = 1:n_reps
            #########################
            rand!(tree, retries=typemax(Int64))
            expr = get_expr(tree)
            f = to_function(problem, expr)
            f(D_1)
            #########################
        end
        t_elapsed = CPUtime_elapsed_s(tstart)
        @notify_observer(observer, "to_function2", [trial, t_elapsed])
    end

    if b_emulated_fitness
        #time emulated_fitness
        println("timing emulated_fitness")
        for trial = 1:n_trials
            tstart = CPUtime_start()
            for i = 1:n_reps
                #########################
                rand!(tree, retries=typemax(Int64))
                expr = get_expr(tree)
                predicts = get_predicts(problem, expr)
                #########################
            end
            t_elapsed = CPUtime_elapsed_s(tstart)
            @notify_observer(observer, "emulated_fitness", [trial, t_elapsed])
        end
    end

    if b_predicts
        #time predicts
        println("timing predicts")
        for trial = 1:n_trials
            rand!(tree, retries=typemax(Int64))
            expr = get_expr(tree)
            println("trial=$trial, expr=", string(expr))
            tstart = CPUtime_start()
            for i = 1:n_reps
                #########################
                predicts = get_predicts(problem, expr)
                #########################
            end
            t_elapsed = CPUtime_elapsed_s(tstart)
            @notify_observer(observer, "predicts", [trial, t_elapsed])
        end
    end

    if b_emulated_predicts
        #time emulated_predicts
        println("timing emulated_predicts")
        for trial = 1:n_trials
            rand!(tree, retries=typemax(Int64))
            expr = get_expr(tree)
            println("trial=$trial, expr=", string(expr))
            tstart = CPUtime_start()
            for i = 1:n_reps
                #########################
                Dl = problem.Dl
                f = to_function(problem, expr)
                predicts = Array(Bool, length(Dl))
                for i = 1:length(Dl)
                    predicts[i] = f(getrecords(Dl, i))
                end
                #########################
            end
            t_elapsed = CPUtime_elapsed_s(tstart)
            @notify_observer(observer, "emulated_predicts", [trial, t_elapsed])
        end
    end

    #time emulated_predicts2
    println("timing emulated_predicts2")
    for trial = 1:n_trials
        rand!(tree, retries=typemax(Int64))
        expr = get_expr(tree)
        Dl = problem.Dl
        f = to_function(problem, expr)
        println("trial=$trial, expr=", string(expr))
        tstart = CPUtime_start()
        for i = 1:n_reps
            #########################
            predicts = Array(Bool, length(Dl))
            for i = 1:length(Dl)
                predicts[i] = f(getrecords(Dl, i))
            end
            #########################
        end
        t_elapsed = CPUtime_elapsed_s(tstart)
        @notify_observer(observer, "emulated_predicts2", [trial, t_elapsed])
    end

    predicts = get_predicts(problem, expr)
    #time metrics
    println("timing metrics")
    for trial = 1:n_trials
        tstart = CPUtime_start()
        for i = 1:n_reps
            #########################
            _, _, metric = gini_metrics(predicts, problem.Dl.labels)
            #########################
        end
        t_elapsed = CPUtime_elapsed_s(tstart)
        @notify_observer(observer, "metrics", [trial, t_elapsed])
    end

    save_log("timingresults.txt", logs)
    logs
end


text:
using StatsBase

"""
the eigenvectors have arbitrary directions.
"""
function test_proj(v1, v2; atol = 0.0)
    all(isapprox(v1[:, i],  v2[:, i], atol = atol) ||
        isapprox(v1[:, i], -v2[:, i], atol = atol) for i in 1:size(v1, 2))
end

@testset "PCA" begin

    x = rand(100, 4)
    w = ones(100)
    zsm = fit(ZScoreTransform, x')
    zm = fit(ZScoreTransform, x', scale = false)

    c = fit!(WeightedCovMatrix(), x, w)

    t, p = pca(c, cov_pca = true)
    p6 = fit(
        PCA,
        StatsBase.transform(zm, permutedims(x)),
        method = :cov,
        maxoutdim = 4,
        pratio = 1.0,
        mean = nothing
    )

    @test t.mean ≈ zm.mean
    @test t.scale ≈ zm.scale
    @test isapprox(p.prinvars, p6.prinvars, atol = 1.0e-2)
    @test isapprox(p.tprinvar, p6.tprinvar, atol = 1.0e-2)
    @test isapprox(p.tvar, p6.tvar, atol = 1.0e-2)
    @test test_proj(p.proj, p6.proj)

    t, p = pca(c, cov_pca = false)
    p4 = fit(
        PCA,
        StatsBase.transform(zsm, x'),
        method = :cov,
        maxoutdim = 4,
        pratio = 1.0,
        mean = 0
    )

    @test t.mean ≈ zsm.mean
    @test isapprox(t.scale, zsm.scale, atol = 1.0e-2)
    @test isapprox(p.prinvars, p4.prinvars, atol = 1.0e-2)
    @test isapprox(p.tprinvar, p4.tprinvar, atol = 1.0e-2)
    @test isapprox(p.tvar, p4.tvar, atol = 1.0e-2)
    @test test_proj(p.proj, p4.proj)

end


text:
@testset "Node Building" begin
    n = Node()
    @test isa(n, AbstractNode)
    @test isa(n, GeneralNode{Float64, Int64})

    n = Node("Name")
    @test isa(n, AbstractNode)
    @test isa(n, GeneralNode{Float64, Int64})
    @test n.name == "Name"

    n = Node("Name", 0.75)
    @test isa(n, AbstractNode)
    @test isa(n, GeneralNode{Float64, Int64})
    @test n.name == "Name"
    @test n.inc_length == 0.75
end

@testset "add_child!" begin
    tree = ParseNewick("(A,B,(C,D,E)F)G;")
    to_add = ParseNewick("no_name;")
    root = find_by_name(tree,"G")
    add_child!(root,to_add)
    @test root.children[end].name == "no_name"
    @test length(post_order(tree))==8

    tree = ParseNewick("(A,B,(C,D,E)F)G;")
    to_add = ParseNewick("no_name;")
    leaf = find_by_name(tree,"D")
    add_child!(leaf,to_add)
    @test leaf.children[end].name == "no_name"
    @test length(post_order(tree))==8

    tree = ParseNewick("(A,B,(C,D,E)F)G;")
    to_add = ParseNewick("no_name;")
    mid = find_by_name(tree,"F")
    add_child!(mid,to_add)

    @test mid.children[end].name == "no_name"
    @test length(post_order(tree))==8
end

@testset "remove_child!" begin
    tree = ParseNewick("((A,B)C,(D,E)F)G;")
    mother = find_by_name(tree,"C")
    remove_child!(mother,true)
    @test length(mother.children) == 1
    @test mother.children[1].name == "B"
    @test length(post_order(tree)) == 6

    tree = ParseNewick("((A,B)C,(D,E)F)G;")
    mother = find_by_name(tree,"C")
    remove_child!(mother,false)
    @test length(mother.children) == 1
    @test mother.children[1].name == "A"
    @test length(post_order(tree)) == 6

    tree = ParseNewick("((A,B)C,(D,E)F)G;")
    remove_child!(tree,true)
    @test tree.children[1].name == "F"
    @test length(tree.children) == 1
    @test length(post_order(tree)) == 4

end

@testset "remove_child! by node" begin
    tree = ParseNewick("((A,B)C,(D,E)F)G;")
    mother = find_by_name(tree,"C")
    to_remove = find_by_name(tree,"A")
    remove_child!(mother,to_remove)
    @test length(mother.children) == 1
    @test mother.children[1].name == "B"
    @test length(post_order(tree)) == 6

    tree = ParseNewick("((A,B)C,(D,E)F)G;")
    mother = find_by_name(tree,"C")
    to_remove = find_by_name(tree,"B")
    remove_child!(mother,to_remove)
    @test length(mother.children) == 1
    @test mother.children[1].name == "A"
    @test length(post_order(tree)) == 6

    tree = ParseNewick("((A,B)C,(D,E)F)G;")
    to_remove = find_by_name(tree,"C")
    remove_child!(tree,to_remove)
    @test tree.children[1].name == "F"
    @test length(tree.children) == 1
    @test length(post_order(tree)) == 4
end

@testset "tree_length" begin
    tree = ParseNewick("((A:5,B:5)C:3.5,(D:5,E:5)F:5)G:5;")
    @test tree_length(tree) == 28.5
end

@testset "tree_height" begin
    tree = ParseNewick("((A:5,B:5)C:3.5,(D:5,E:5)F:5)G:5;")
    @test tree_height(tree) == 10.0
end

@testset "node_height" begin
    tree = ParseNewick("((A:5,B:5)C:9,(D:5,E:5)F:5)G:5;")
    node_height(tree)
    @test tree.height == 14.0
end

@testset "node_age" begin
    tree = ParseNewick("(((A:8,B:5)F:2,C:10)G:1,(D:12,E:4)H:3)R:1;")
    leaf = find_by_name(tree, "A")
    internal_node = find_by_name(tree, "F")
    furthest_leaf = find_by_name(tree, "D")
    @test node_age(furthest_leaf) == 0.0
    @test node_age(leaf) == 4.0
    @test node_age(internal_node) == 12.0
    @test node_age(tree) == 15.0
end

@testset "node_depth" begin
    tree = ParseNewick("((A:5,B:5)C:9,(D:5,E:5)F:5)G:5;")
    mid = find_by_name(tree,"F")
    leaf = find_by_name(tree, "A")
    @test node_depth(tree) == 0
    @test node_depth(mid) == 1
    @test node_depth(leaf) == 2
end

@testset "get_path" begin
    tree = ParseNewick("((A:5,B:5)C:9,(D:5,E:5)F:5)G:5;")
    mid = find_by_name(tree,"F")
    leaf = find_by_name(tree,"A")
    @test get_path(tree,mid) == [6]
    @test get_path(tree,leaf) == [1, 5]
end

@testset "path_length" begin
    tree = ParseNewick("((A:5,B:5)C:9,(D:5,E:5)F:5)G:5;")
    mid = find_by_name(tree,"F")
    leaf = find_by_name(tree,"A")
    @test path_length(tree,mid) == 5.0
    @test path_length(tree,leaf) == 14.0
end

@testset "get_sister" begin
    tree = ParseNewick("((A:5,B:5)C:9,(D:5,E:5)F:5)G:5;")
    mid = find_by_name(tree,"F")
    leaf = find_by_name(tree,"A")
    leafresult = get_sister(leaf)
    @test leafresult.name == "B"
    midresult = get_sister(mid)
    @test midresult.name == "C"
end

@testset "set_binary!" begin
    tree = parsing_newick_string("((A:5,B:5)C:9,(D:5,E:5)F:5)G:5;")
    MCPhyloTree.set_binary!(tree)
    mid = find_by_name(tree,"F")
    leaf = find_by_name(tree,"A")
    @test tree.binary == "1"
    @test leaf.binary =="1,0,0"
    @test mid.binary == "1,1"
end

@testset "number_nodes!" begin
    tree = parsing_newick_string("((B:5,A:5)C:9,(D:5,E:5)F:5)G:5;")
    MCPhyloTree.number_nodes!(tree)
    mid = find_by_name(tree,"F")
    leaf = find_by_name(tree,"A")
    @test tree.num == 7
    @test mid.num == 6
    @test leaf.num == 1
end

@testset "blv" begin
    tree = ParseNewick("((B:5,A:5)C:9,(D:5,E:5)F:5)G:5;")
    @test get_branchlength_vector(tree) == [5.0, 5.0, 5.0, 5.0, 9.0, 5.0]
    set_branchlength_vector!(tree,[1.0,1.0,1.0,1.0,1.0,1.0])
    @test get_branchlength_vector(tree) == [1.0,1.0,1.0,1.0,1.0,1.0]
    @test MCPhyloTree.get_sum_seperate_length!(tree) == [2.0, 4.0, 0.0, 0.0]
end

@testset "internal_external" begin
    tree = ParseNewick("((B:5,A:5)C:9,(D:5,E:5)F:5)G:5;")
    map = MCPhyloTree.internal_external_map(tree)
    @test map == [0, 0, 0, 0, 1, 1]
    @test MCPhyloTree.internal_external(tree) == [0, 0, 0, 0, 1, 1]
end

@testset "check_binary" begin
    tree = ParseNewick("((B:5,A:5)C:9,(D:5,E:5)F:5)G:5;")
    @test check_binary(tree)
    tree = ParseNewick("((B:5)C:9,(D:5,E:5)F:5)G:5;")
    @test !check_binary(tree)
    tree = ParseNewick("((B:5,A:5,X:5)C:9,(D:5,E:5)F:5)G:5;")
    @test !check_binary(tree)
end

@testset "insert_node!" begin
    tree = ParseNewick("(A,B,(C,D,E)F)G;")
    insert_tree = ParseNewick("(A,B,((C,D,E)no_name)F)G;")
    insert_tree_newick = newick(insert_tree)
    insert_tree2 = ParseNewick("(A,B,(((C,D)no_name,E)no_name)F)G;")
    insert_tree_newick2 = newick(insert_tree2)

    children = [find_by_name(tree, "C"), find_by_name(tree, "D"),
                find_by_name(tree, "E")]
    mother = find_by_name(tree, "F")
    inserted_node = insert_node!(mother, children)
    @test newick(tree) == insert_tree_newick

    pop!(children)
    insert_node!(inserted_node, children)
    @test newick(tree) == insert_tree_newick2

    new_node = Node("G")
    append!(children, new_node)
    @test_throws AssertionError insert_node!(mother, children)
 end

@testset "delete_node!" begin
    tree = ParseNewick("(A,B,(C,D)E)F;")
    remove_tree = ParseNewick("(A,B,C,D)F;")
    remove_tree_newick = newick(remove_tree)
    remove_tree2 = ParseNewick("(B,C,D)F;")
    remove_tree2_newick = newick(remove_tree2)

    delete_node!(find_by_name(tree, "E"))
    @test newick(tree) == remove_tree_newick
    delete_node!(find_by_name(tree, "A"))
    @test newick(tree) == remove_tree2_newick
    @test_throws ArgumentError delete_node!(find_by_name(tree, "F"))
end

@testset "find_lca" begin
    tree = ParseNewick("(A,B,(C,(D,E)F)G)H;")

    @testset "find_lca by name" begin
        @test find_lca(tree, ["A","C"]) == tree
        @test find_lca(tree, ["D", "E"]) == find_by_name(tree, "F")
        @test find_lca(tree, ["D", "F"]) == find_by_name(tree, "F")
        @test find_lca(tree, ["C", "D"]) == find_by_name(tree, "G")
        @test find_lca(tree, ["A","B","C"]) == tree
        @test find_lca(tree, ["D", "E", "F", "G"]) == find_by_name(tree, "G")
        @test find_lca(tree, ["D", "E", "C"]) == find_by_name(tree, "G")
    end # find_lca by name

    @testset "find_lca by nodes" begin
        nodes = [find_by_name(tree, "A"), find_by_name(tree, "B"),
                 find_by_name(tree, "C")]
        @test find_lca(tree, nodes) == tree

        nodes = [find_by_name(tree, "D"), find_by_name(tree, "E"),
                 find_by_name(tree, "F"), find_by_name(tree, "G")]
        @test find_lca(tree, nodes) == find_by_name(tree, "G")
    end # find_lca by nodes
end

@testset "check_leafsets" begin
    tree1 = ParseNewick("(((A,B),C),(D,E));")
    tree2 = ParseNewick("((A,C),(B,D,E));")
    tree3 = ParseNewick("(((G,C),A),D,F);")
    tree4 = ParseNewick("(((B,C),A),D,E);")
    tree5 = ParseNewick("(((B,C),A),D,F);")

    trees = [tree1, tree2, tree3, tree4, tree5]

    @test_throws ArgumentError check_leafsets(trees)
end

@testset "mother" begin
    tree1 = ParseNewick("(((A,B)F,C)G,(D,E)H)R;")
    @test_throws ArgumentError get_mother(find_by_name(tree1, "R"))
    @test get_mother(find_by_name(tree1, "A")).name == "F"
end

@testset "force_ultrametric" begin
    tree1 = ParseNewick("(((A:5,B:5)F:5,C:5)G:5,(D:5,E:5)H:5)R:5;")
    MCPhyloTree.force_ultrametric!(tree1)
    for x in get_leaves(tree1)
        @test path_length(tree1,x) == 15.0
    end
end

@testset "ascii" begin
    tree = ParseNewick("(A,B)C;")
    singletree = ParseNewick("(((A,B)C)D,(E)F,G)H;")
    bigtree = ParseNewick("((A,B)C,(D,E)F)G;")
    lines,_ = MCPhyloTree.ascii(tree)
    biglines,_ = MCPhyloTree.ascii(bigtree)
    singlelines,_ = MCPhyloTree.ascii(singletree)
    @test lines == ["   /-A", "-C|", "   \\-B"]
    @test biglines == ["      /-A", "   /C|", "  |   \\-B", "-G|", "  |   /-D", "   \\F|", "      \\-E"]
    @test singlelines == ["         /-A", "   /D -C|", "  |      \\-B", "-H|", "  |-F --E", "  |", "   \\-G"]
end


text:
include("HWLC.jl")
using .HWLC
using Plots, Distributions, MAT

#Set plot parameters
gr(xlabel="F(X) (pN)", label="",ylabel="F'(X) (pN/nm)",
    grid=false,guidefont=("Arial",10),
    legend=:topleft,
    legendfont=("Arial",10),scale=:log10)

#Import experiment data for comparison
data=matread("Experimental_data/U2OS_ForceStiffnessStd.mat")
mean_experiment=vec(data["k_avg"])
std_exp=vec(data["k_err"])
fs_exp=vec(data["f_ax"])

#Set parameters
kT=4.22
Marko_Siggia=true

f_min,f_max,N_f=extrema(fs_exp)...,500
forces=exp.(range(log(f_min),log(f_max),length=N_f))

N=30        #Number of component chains
L=3000.     #Total length of assembly

#Power-law distribution for persistence lengths
slope=0.86        #Slope of stiffness-force graph
α= slope -2.     #Power law P(L_p)~L_p^(α)
minp, maxp=kT/100000,kT/9.

"""
Show individual and total curves for one assembly
"""
persistence_lengths=power_law_samples(α,N,minp,maxp)
#ls=ones(N)*L/N  #All equal length
ls=rand(Normal(L/N,L/N*0.1),N) #Alternatively sample from normal distr

total_ext,extensions=force_extensions(forces,ls,persistence_lengths,kT,Marko_Siggia)
total_stiffness=stiffnesses(forces,total_ext)[1]
individual_stiffnesses=stiffnesses(forces,extensions)

#Plot individual and total response
plot(forces[2:end],total_stiffness, label="Assembly");
plot!([20,200],[10^-1.5,(200/20)^slope*10^-1.5], color="black",label="Slope $slope");
for ind in individual_stiffnesses
    plot!(forces[2:end],ind, color="lightgrey")
end
display(current())
png("Plots/power_law_individual_responses")

"""
Sample many assemblies, and compare mean
stiffening behaviour to data
"""
n_samples=15
sim_data=zeros(length(forces)-1,n_samples)
for i in 1:n_samples
    #sample assembly
    sample_ls=rand(Normal(L/N,L/N*0.1),N) #comment if using constant lengths
    sample_ps=power_law_samples(α,N,minp,maxp)
    tot,exts=force_extensions(forces,sample_ls,sample_ps,kT,Marko_Siggia)
    sim_data[:,i]=stiffnesses(forces,tot)[1]
end

mean_sim=mean(sim_data,dims=2)
std_sim=std(sim_data,dims=2)

plot(fs_exp,mean_experiment, label="Experimental mean",legend=:topleft, color="blue");
plot!(fs_exp,mean_experiment.+std_exp, color="lightblue");
plot!(fs_exp[findall(>(0.),mean_experiment.-std_exp)],(mean_experiment.-std_exp)[findall(>(0.),mean_experiment.-std_exp)], color="lightblue");
plot!(forces[2:end],mean_sim, label="Sim mean", color="orange", ribbon=std_sim)
png("Plots/power_law_compare_data_sim")


text:
include("src\\ConceptDrift.jl")

a = ConceptDrift.SineGenerator(1)
c = ConceptDrift.NaiveBayes()
for i in 1:1000000
    X,y = ConceptDrift.next_sample!(a, 1)
    ConceptDrift.partial_fit!(c, vec(X), y[1])
end

right = 0
wrong = 0
for i in 1:100000
    X,y = ConceptDrift.next_sample!(a, 1)
    label = ConceptDrift.predict(c, vec(X))
    if label == y[1]
        right += 1
    else
        wrong += 1
    end
end

acc = right / (right + wrong)


text:
#=

  From BLOG example/burglary-multihouse.blog

  Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
       Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

  burglary[1]    0.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
  burglary[2]    0.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
  burglary[3]    0.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
  burglary[4]    0.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
   earthquake    0.9208    0.2701     0.0027    0.0270     23.2091    1.0899        0.5575
     alarm[1]    0.9994    0.0245     0.0002    0.0006   1543.6492    1.0005       37.0775
     alarm[2]    0.9992    0.0283     0.0003    0.0008   1153.7201    1.0007       27.7117
     alarm[3]    0.9992    0.0283     0.0003    0.0008   1153.7201    1.0007       27.7117
     alarm[4]    0.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN

  Cf ~/webppl/burglary_multihouse.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function burglary_multihouse(evidence=[true,true,true,false])
    mary = 1
    john = 2
    cathy = 3
    roger = 4
    n = 4

    burglary ~ filldist(flip(0.003),n)
    earthquake ~ flip(0.002)
    alarm = tzeros(n)
    for p in 1:n
        alarm[p] ~ (burglary[p] == false && earthquake ==false)  ? flip(0.01) :
                   (burglary[p] == false && earthquake == true)  ? flip(0.40) :
                   (burglary[p] == true  && earthquake == false) ? flip(0.80) :
                   (burglary[p] == true  && earthquake == true ) ? flip(0.90) : flip(0.0)
    end

    # for p in 1:n
    #     true ~ Dirac(alarm[p] == evidence[p])
    # end
    true ~ Dirac(alarm[mary] == true)
    true ~ Dirac(alarm[john] == true)
    true ~ Dirac(alarm[cathy] == true)
    true ~ Dirac(alarm[roger] == false)
end

# alarms from [mary,john,cathy,roger]
evidence=[true,true,true,false]
model = burglary_multihouse(evidence)

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(10_000), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)
# chns = sample(model, HMC(0.1,10), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#  Knitro example defining a single callback function for function
#  and gradient(i.e. EVALFCGA) requests, and a second callback for
#  Hessian(i.e. EVALH / EVALHV) requests.  This model is test
#  problem HS40 from the Hock & Schittkowski collection.
#
#  max   x0*x1*x2*x3        (obj)
#  s.t.  x0^3 + x1^2 = 1    (c0)
#        x0^2*x3 - x2 = 0   (c1)
#        x3^2 - x1 = 0      (c2)
#
#  See exampleNLP2.c to see how this is solved in the standard
#  way with separate callbacks for functions and gradients.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO
using Test

function example_fcga(; verbose=true)
    #*------------------------------------------------------------------*/
    #*     FUNCTION callbackEvalFCGA                                    */
    #*------------------------------------------------------------------*/
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # To compute the functions and gradient together, set "obj", "c",
    # "objGrad" and "jac" in the KNITRO.KN_eval_result structure when there
    # is a function+gradient evaluation request(i.e. EVALFCGA).
    #
    # NOTE: It is generally more efficient and recommended to have
    # separate callback routines for functions and gradients since a
    # gradient evaluation is not always needed for every function
    # evaluation.  However, in some cases, it may be more convenient
    # to compute them together if most of the work for computing the
    # gradients is already done for the function evaluation.
    function callbackEvalFCGA(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate nonlinear term in objective
        evalResult.obj[1] = x[1] * x[2] * x[3] * x[4]

        # Evaluate nonlinear terms in constraints
        evalResult.c[1] = x[1] * x[1] * x[1]
        evalResult.c[2] = x[1] * x[1] * x[4]

        # Evaluate nonlinear term in objective gradient
        evalResult.objGrad[1] = x[2] * x[3] * x[4]
        evalResult.objGrad[2] = x[1] * x[3] * x[4]
        evalResult.objGrad[3] = x[1] * x[2] * x[4]
        evalResult.objGrad[4] = x[1] * x[2] * x[3]

        # Evaluate nonlinear terms in constraint gradients(Jacobian)
        evalResult.jac[1] = 3.0 * x[1] * x[1] # derivative of x0^3 term  wrt x0
        evalResult.jac[2] = 2.0 * x[1] * x[4] # derivative of x0^2*x3 term  wrt x0
        evalResult.jac[3] = x[1] * x[1]     # derivative of x0^2*x3 terms wrt x3

        return 0
    end

    #*------------------------------------------------------------------*/
    #*     FUNCTION callbackEvalH                                       */
    #*------------------------------------------------------------------*/
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "hess" and "hessVec" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalH(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x
        lambda_ = evalRequest.lambda
        # Scale objective component of Hessian by sigma
        sigma = evalRequest.sigma

        # Evaluate nonlinear terms in the Hessian of the Lagrangian.
        # Note: If sigma=0, some computations can be avoided.
        if sigma > 0.0  # Evaluate the full Hessian of the Lagrangian
            evalResult.hess[1] = lambda_[1] * 6.0 * x[1] + lambda_[2] * 2.0 * x[4]
            evalResult.hess[2] = sigma * x[3] * x[4]
            evalResult.hess[3] = sigma * x[2] * x[4]
            evalResult.hess[4] = sigma * x[2] * x[3] + lambda_[2] * 2.0 * x[1]
            evalResult.hess[5] = sigma * x[1] * x[4]
            evalResult.hess[6] = sigma * x[1] * x[3]
            evalResult.hess[7] = sigma * x[1] * x[2]
        else            # sigma=0, do not include objective component
            evalResult.hess[1] = lambda_[1] * 6.0 * x[1] + lambda_[2] * 2.0 * x[4]
            evalResult.hess[2] = 0.0
            evalResult.hess[3] = 0.0
            evalResult.hess[4] = lambda_[2] * 2.0 * x[1]
            evalResult.hess[5] = 0.0
            evalResult.hess[6] = 0.0
            evalResult.hess[7] = 0.0
        end

        return 0
    end

    #*------------------------------------------------------------------*
    #*     main                                                         *
    #*------------------------------------------------------------------*

    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Initialize Knitro with the problem definition.

    # Add the variables and specify initial values for them.
    # Note: any unset lower bounds are assumed to be
    # unbounded below and any unset upper bounds are
    # assumed to be unbounded above.
    vars = KNITRO.KN_add_vars(kc, 4)
    for x in vars
        KNITRO.KN_set_var_primal_init_value(kc, x, 0.8)
    end

    # Add the constraints and set the rhs and coefficients
    KNITRO.KN_add_cons(kc, 3)
    KNITRO.KN_set_con_eqbnds_all(kc, [1.0, 0.0, 0.0])

    # Coefficients for 2 linear terms
    lconIndexCons = Int32[1, 2]
    lconIndexVars = Int32[2, 1]
    lconCoefs = [-1.0, -1.0]
    KNITRO.KN_add_con_linear_struct(kc, lconIndexCons, lconIndexVars, lconCoefs)

    # Coefficients for 2 quadratic terms

    qconIndexCons = Int32[0, 2]
    qconIndexVars1 = Int32[1, 3]
    qconIndexVars2 = Int32[1, 3]
    qconCoefs = [1.0, 1.0]

    KNITRO.KN_add_con_quadratic_struct(
        kc,
        qconIndexCons,
        qconIndexVars1,
        qconIndexVars2,
        qconCoefs,
    )

    # Add callback to evaluate nonlinear(non-quadratic) terms in the model:
    #    x0*x1*x2*x3  in the objective
    #    x0^3         in first constraint c0
    #    x0^2*x3      in second constraint c1
    cIndices = Int32[0, 1]
    cb = KNITRO.KN_add_eval_callback(kc, true, cIndices, callbackEvalFCGA)

    # Set obj. gradient and nonlinear jac provided through callbacks.
    # Mark objective gradient as dense, and provide non-zero sparsity
    # structure for constraint Jacobian terms. Set the gradient
    # evaluations to use the same callback routine as used for
    # function evaluations.
    cbjacIndexCons = Int32[0, 1, 1]
    cbjacIndexVars = Int32[0, 0, 3]
    KNITRO.KN_set_cb_grad(
        kc,
        cb,
        nothing,
        jacIndexCons=cbjacIndexCons,
        jacIndexVars=cbjacIndexVars,
    )

    # Set nonlinear Hessian provided through callbacks. Since the
    # Hessian is symmetric, only the upper triangle is provided.
    # The upper triangular Hessian for nonlinear callback structure is:
    #
    #  lambda0*6*x0 + lambda1*2*x3     x2*x3    x1*x3    x1*x2 + lambda1*2*x0
    #               0                    0      x0*x3         x0*x2
    #                                             0           x0*x1
    #                                                          0
    # (7 nonzero elements)
    cbhessIndexVars1 = Int32[0, 0, 0, 0, 1, 1, 2]
    cbhessIndexVars2 = Int32[0, 1, 2, 3, 2, 3, 3]
    KNITRO.KN_set_cb_hess(
        kc,
        cb,
        length(cbhessIndexVars1),
        callbackEvalH,
        hessIndexVars1=cbhessIndexVars1,
        hessIndexVars2=cbhessIndexVars2,
    )

    # Set minimize or maximize(if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MAXIMIZE)

    # Set option to print output after every iteration.
    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ITER : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Set option to tell Knitro that the gradients are being provided
    # with the functions in one callback.
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_EVAL_FCGA, KNITRO.KN_EVAL_FCGA_YES)

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)

    # An example of obtaining solution information.
    if verbose
        println("Knitro converged with final status = ", nStatus)
        println("  optimal objective value  = ", objSol)
        println("  optimal primal values x  = ", x)
        println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
        println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))
    end

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example HS40 FCGA" begin
        @test nStatus == 0
        @test objSol ≈ 0.25
        @test x ≈ [0.793701, 0.707107, 0.529732, 0.840896] atol = 1e-5
    end
end

example_fcga(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_volumerepresentation

"""
    vtk_volumerepresentation(;kwargs...)
    vtk_volumerepresentation(children::Any;kwargs...)
    vtk_volumerepresentation(children_maker::Function;kwargs...)

A VolumeRepresentation component.
VolumeRepresentation is responsible to convert a vtkPolyData into rendering
It takes the following set of properties:
  - colorBy: ['POINTS', ''],
  - pointSize: 1,
  - color: [1,1,1],
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component.
- `colorDataRange` (Array of Reals | String; optional): Data range use for the colorMap
- `colorMapPreset` (String; optional): Preset name for the lookup table color map
- `mapper` (Dict; optional): Properties to set to the mapper
- `property` (Dict; optional): Properties to set to the volume.property
- `volume` (Dict; optional): Properties to set to the volume
"""
function vtk_volumerepresentation(; kwargs...)
        available_props = Symbol[:children, :id, :colorDataRange, :colorMapPreset, :mapper, :property, :volume]
        wild_props = Symbol[]
        return Component("vtk_volumerepresentation", "VolumeRepresentation", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_volumerepresentation(children::Any; kwargs...) = vtk_volumerepresentation(;kwargs..., children = children)
vtk_volumerepresentation(children_maker::Function; kwargs...) = vtk_volumerepresentation(children_maker(); kwargs...)

max_stars_repo_name: rodrigomologni/dash-vtk

max_issues_repo_name: rodrigomologni/dash-vtk

max_forks_repo_name: rodrigomologni/dash-vtk


text:
include("../src/helper.jl")
using LaTeXStrings
using Plots

function circle(h,k,r)
    th = LinRange(0,2*pi,500)
    h .+ r*sin.(th), k.+ r*cos.(th)
end

function plot_hypergraph(xy,Edges,radius,fsize,lw,al,gverts,bverts)

    k = size(Edges,2)
    n = size(xy,1)
    xy_tuples = Vector{Tuple{Float64,Float64}}()
    for i = 1:size(xy,1)
        push!(xy_tuples,(xy[i,1],xy[i,2]))
    end

    # Hyperedges as triangles
    p = plot(grid = false,legend = false,axis = false)
    for i = 1:size(Edges,1)
        sp = Shape(xy_tuples[Edges[i,:]])
        plot!(sp,color = :black, alpha = al,linewidth = lw)
    end

    # Now vertices
    for i = gverts
    plot!(p,circle(xy[i,1],xy[i,2],radius),linewidth = 0,aspect_ratio = 1,fillalpha = 1, c = :green,seriestype = [:shape])
    annotate!(p,xy[i,1],xy[i,2],text("$i",fsize,:white))
    end

    for i = bverts
    plot!(p,circle(xy[i,1],xy[i,2],radius),linewidth = 0,aspect_ratio = 1,fillalpha = 1, c = :blue,seriestype = [:shape])
    annotate!(p,xy[i,1],xy[i,2],text("$i",fsize,:white))
    end

    return p

end

function plot_bipartite(a,dg,db,dif,Edges,radius,grad,fsize,lw,gverts,bverts)
p1 = plot(grid = false,legend = false,axis = false)

xy =[-a 3*dg;
-a dg;
-a -dg;
-a -3*dg;
a db;
a 0;
a -db
]

m = size(Edges,1)
spn = m*dif/2
xyg =[ zeros(m,1) LinRange(spn,-spn,m)]

# Plot lines for edges
for j = 1:m
    edge = Edges[j,:]
    gp_xy = xyg[j,:]
    for i = 1:3
        nd = edge[i]
        if nd > 4
            cl = :blue
        else
            cl = :green
        end
        nd_xy = xy[edge[i],:]
        plot!(p1,[gp_xy[1];nd_xy[1]],[gp_xy[2];nd_xy[2]], linewidth = lw, linecolor = cl)
    end
end

# Nodes on top
n = size(xy,1)
alphabet = ["a";"b";"c";"d";"e";"f";"g";"h";"i";"j";"k";"l"]
# Plot
xyg_tuples = Vector{Tuple{Float64,Float64}}()
for i = 1:size(xyg,1)
    push!(xyg_tuples,(xyg[i,1],xyg[i,2]))
end
xy_tuples = Vector{Tuple{Float64,Float64}}()
for i = 1:size(xy,1)
    push!(xy_tuples,(xy[i,1],xy[i,2]))
end

# Now vertices
for i = 1:n
    if in(i,gverts)
        plot!(p1,circle(xy[i,1],xy[i,2],radius),linewidth = 0,aspect_ratio = 1,fillalpha = 1, c = :green,seriestype = [:shape])
        annotate!(p1,xy[i,1],xy[i,2],text("$i",fsize,:white))
    else
        plot!(p1,circle(xy[i,1],xy[i,2],radius),linewidth = 0,aspect_ratio = 1,fillalpha = 1, c = :blue,seriestype = [:shape])
        annotate!(p1,xy[i,1],xy[i,2],text("$i",fsize,:white))
    end
end

# Group Vertices
for i = 1:m
plot!(p1,circle(xyg[i,1],xyg[i,2],grad),linewidth = 0,aspect_ratio = 1,fillalpha = 1, c = :gray,seriestype = [:shape])
annotate!(p1,xyg[i,1],xyg[i,2],text("$(alphabet[i])",fsize,:white))
end

return p1
end

function vert_table()
    ## Vertical table
    for i = 1:n
        if i < 5
            print("{\\color{dgreen}$i} & ")
        else
            print("{\\color{blue}$i} & ")
        end
        for j = 1:2
            print("$(d[i,j]) & ")
        end
        print("$(d[n,3]) \\\\\n")
    end
end

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
function fit(::Type{LogNormal}, m::AbstractMoments)
    # https://en.wikipedia.org/wiki/Log-normal_distribution
    n_moments(m) >= 2 || error("Need mean and variance to estimate lognormal")
    γ = 1+var(m)/mean(m)^2
    μ = log(mean(m)/sqrt(γ))
    σ = sqrt(log(γ))
    return LogNormal(μ,σ)
end

function fit(::Type{LogNormal}, lower::QuantilePoint, upper::QuantilePoint)
    #length(qset) == 2 || error("only implemented yet for exactly two quantiles.")
    #qset_log = [QuantilePoint(qp, q = log(qp.q)) for qp in qset]
    lower_log = QuantilePoint(lower, q = log(lower.q))
    upper_log = QuantilePoint(upper, q = log(upper.q))
    DN = fit(Normal, lower_log, upper_log)
    LogNormal(params(DN)...)
end

function fit_mean_quantile(::Type{LogNormal}, mean::Real, qp::QuantilePoint)
    # solution of
    # (1) mean = exp(mu + sigma^2/2)
    # (2) upper = mu + sigmaFac sigma
    # see R packaage lognorm inst/doc/coefLognorm.Rmd for derivation
    sigmaFac = quantile(Normal(),qp.p)
    m = log(mean)
    discr = sigmaFac^2 - 2*(log(qp.q) - m)
    (discr > 0) || error("Cannot fit LogNormal with quantile $(qp) and mean $(mean).")
    sigma = sigmaFac > sqrt(discr) ? (sigmaFac - sqrt(discr)) : sigmaFac + sqrt(discr)
    mu = m - sigma^2/2
    LogNormal(mu, sigma)
end

function fit_mode_quantile(::Type{LogNormal}, mode::Real, qp::QuantilePoint)
    # solution of
    # (1) mle = exp(mu - sigma^2)
    # (2) upper = mu + sigmaFac sigma
    # see R packaage lognorm inst/doc/coefLognorm.Rmd for derivation
    sigmaFac = quantile(Normal(),qp.p)
    m = log(mode)
    discr = sigmaFac^2/4 + (log(qp.q) - m)
    (discr > 0) || error("Cannot fit LogNormal with quantile $(qp) and mode $(mode).")
    root_discr = sqrt(discr)
    mfh = -sigmaFac/2
    sigma = mfh > root_discr ? (mfh - root_discr) : (mfh + root_discr)
    #sigma = mfh + root_discr
    mu = m + sigma^2
    LogNormal(mu, sigma)
end

"""
    Σstar <: AbstractΣstar

Represent the multiplicative standard deviation of a LogNormal distribution.

Supports dispatch of `fit`.
Invoking the type as a function returns its single value.

# Examples
```jldoctest; output = false, setup = :(using DistributionFits)
a = Σstar(4.2)
a() == 4.2
# output
true
```
"""
abstract type AbstractΣstar end

struct Σstar <: AbstractΣstar
    σstar
end
(a::Σstar)() = a.σstar

"""
    σstar(d)

Get the multiplicative standard deviation of LogNormal distribution d.

# Arguments
- `d`: The type of distribution to fit

# Examples
```jldoctest fm1; output = false, setup = :(using DistributionFits)
d = LogNormal(2,log(1.2))
σstar(d) == 1.2
# output
true
```
"""
σstar(d::LogNormal) = exp(params(d)[2])

"""
    fit(D, mean, σstar)

Fit a statistical distribution of type `D` to mean and multiplicative
standard deviation.

# Arguments
- `D`: The type of distribution to fit
- `mean`: The moments of the distribution
- `σstar::AbstractΣstar`: The multiplicative standard deviation

See also [`σstar`](@ref), [`AbstractΣstar`](@ref).

# Examples
```jldoctest fm1; output = false, setup = :(using DistributionFits)
d = fit(LogNormal, 2, Σstar(1.1));
(mean(d), σstar(d)) == (2, 1.1)
# output
true
```
"""
function fit(::Type{LogNormal}, mean, σstar::AbstractΣstar)
    σ = log(σstar())
    μ = log(mean) - σ*σ/2
    LogNormal(μ, σ)
end

"""
    fit_mean_relerror(D, mean, relerror)

Fit a distribution of type `D` to mean and relative error.

# Arguments
- `D`: The type of distribution to fit
- `mean`: The first moment of the distribution
- `relerror`: The relative error, i.e. the coefficient of variation

# Examples
```jldoctest fm1; output = false, setup = :(using DistributionFits)
d = fit_mean_relerror(LogNormal, 10.0, 0.03);
(mean(d), std(d)/mean(d)) .≈ (10.0, 0.03)
# output
(true, true)
```
"""
function fit_mean_relerror(::Type{LogNormal}, mean, relerror)
    # e.g. Limpert 2001, Wutzler 2020
    w = 1 + abs2(relerror)
    μ = log(mean/sqrt(w))
    σ = sqrt(log(w))
    LogNormal(μ, σ)
end


text:
trtrs!(c1::Char, c2::Char, c3::Char, r::AbstractMatrix, b::AbstractVecOrMat) = LinearAlgebra.LAPACK.trtrs!(c1, c2, c3, r, b)

"""
    copyltu!(A::AbstractMatrix) -> AbstractMatrix
copy the lower triangular to upper triangular.
"""
function copyltu!(A::AbstractMatrix)
    m, n = size(A)
    for i=1:m
        A[i,i] = real(A[i,i])
        for j=i+1:n
            @inbounds A[i,j] = conj(A[j,i])
        end
    end
    A
end

do_adjoint(A::Matrix) = Matrix(A')


text:
"""
```
plot_shock_decomposition(m, var, class, input_type, cond_type;
    title = "", kwargs...)

plot_shock_decomposition(m, vars, class, input_type, cond_type;
    forecast_string = "", groups = shock_groupings(m),
    plotroot::String = figurespath(m, \"forecast\"), titles = [],
    kwargs...)
```

Plot shock decomposition(s) for `var` or `vars`.

### Inputs

- `m::AbstractModel`
- `var::Symbol` or `vars::Vector{Symbol}`: variable(s) whose shock decomposition
  is to be plotted, e.g. `:obs_gdp` or `[:obs_gdp, :obs_nominalrate]`
- `class::Symbol`
- `input_type::Symbol`
- `cond_type::Symbol`

### Keyword Arguments

- `forecast_string::String`
- `groups::Vector{ShockGroup}`
- `plotroot::String`: if nonempty, plots will be saved in that directory
- `title::String` or `titles::Vector{String}`
- `verbose::Symbol`

See `?shockdec` for additional keyword arguments, all of which can be passed
into `plot_history_and_forecast`.

### Output

- `p::Plot` or `plots::OrderedDict{Symbol, Plot}`
"""
function plot_shock_decomposition(m::AbstractModel, var::Symbol, class::Symbol,
                                  input_type::Symbol, cond_type::Symbol;
                                  title = "",
                                  kwargs...)

    plots = plot_shock_decomposition(m, [var], class, input_type, cond_type;
                                     titles = isempty(title) ? String[] : [title],
                                     kwargs...)
    return plots[var]
end

function plot_shock_decomposition(m::AbstractModel, vars::Vector{Symbol}, class::Symbol,
                                  input_type::Symbol, cond_type::Symbol;
                                  forecast_string::String = "",
                                  groups::Vector{ShockGroup} = shock_groupings(m),
                                  plotroot::String = figurespath(m, "forecast"),
                                  titles::Vector{String} = String[],
                                  verbose::Symbol = :low,
                                  kwargs...)
    # Read in MeansBands
    output_vars = [Symbol(prod, class) for prod in [:shockdec, :trend, :dettrend, :hist, :forecast]]
    mbs = map(output_var -> read_mb(m, input_type, cond_type, output_var, forecast_string = forecast_string),
              output_vars)

    # Get titles if not provided
    if isempty(titles)
        detexify_title = typeof(Plots.backend()) == Plots.GRBackend
        titles = map(var -> describe_series(m, var, class, detexify = detexify_title), vars)
    end

    # Loop through variables
    plots = OrderedDict{Symbol, Plots.Plot}()
    for (var, title) in zip(vars, titles)
        # Call recipe
        plots[var] = shockdec(var, mbs..., groups;
                              ylabel = series_ylabel(m, var, class) * "\n(deviations from mean)",
                              title = title, kwargs...)

        # Save plot
        if !isempty(plotroot)
            output_file = get_forecast_filename(plotroot, filestring_base(m), input_type, cond_type,
                                                Symbol("shockdec_", detexify(var)),
                                                forecast_string = forecast_string,
                                                fileformat = plot_extension())
            save_plot(plots[var], output_file, verbose = verbose)
        end
    end
    return plots
end

@userplot Shockdec

"""
```
shockdec(var, shockdec, trend, dettrend, hist, forecast, groups;
    start_date = shockdec.means[1, :date],
    end_date = shockdec.means[end, :date],
    hist_label = \"Detrended History\",
    forecast_label = \"Detrended Forecast\",
    hist_color = :black, forecast_color = :red, tick_size = 5)
```

User recipe called by `plot_shock_decomposition`.

### Inputs

- `var::Symbol`: e.g. `:obs_gdp`
- `shockdec::MeansBands`
- `trend::MeansBands`
- `dettrend::MeansBands`
- `hist::MeansBands`
- `forecast::MeansBands`
- `groups::Vector{ShockGroup}`

### Keyword Arguments

- `start_date::Date`
- `end_date::Date`
- `hist_label::String`
- `forecast_label::String`
- `hist_color`
- `forecast_color`
- `tick_size::Int`: x-axis (time) tick size in units of years

Additionally, all Plots attributes (see docs.juliaplots.org/latest/attributes)
are supported as keyword arguments.
"""
shockdec

@recipe function f(sd::Shockdec;
                   start_date = sd.args[2].means[1, :date],
                   end_date = sd.args[2].means[end, :date],
                   hist_label = "Detrended History",
                   forecast_label = "Detrended Forecast",
                   hist_color = :black,
                   forecast_color = :red,
                   tick_size = 5)
    # Error checking
    if length(sd.args) != 7 || typeof(sd.args[1]) != Symbol ||
        typeof(sd.args[2]) != MeansBands || typeof(sd.args[3]) != MeansBands ||
        typeof(sd.args[4]) != MeansBands || typeof(sd.args[5]) != MeansBands ||
        typeof(sd.args[6]) != MeansBands || typeof(sd.args[7]) != Vector{ShockGroup}

        error("shockdec must be given a Symbol, five MeansBands, and a Vector{ShockGroup}. Got $(typeof(sd.args))")
    end

    var, shockdec, trend, dettrend, hist, forecast, groups = sd.args

    # Construct DataFrame with detrended mean, deterministic trend, and all shocks
    df = DSGE.prepare_means_table_shockdec(shockdec, trend, dettrend, var,
                                      mb_hist = hist, mb_forecast = forecast,
                                      detexify_shocks = false,
                                      groups = groups)
    dates = df[:date]
    xnums = (1:length(dates)) - 0.5

    # Assign date ticks
    inds = intersect(find(x -> start_date .<= x .<= end_date,  dates),
                     find(x -> Dates.month(x) == 3,            dates),
                     find(x -> Dates.year(x) % tick_size == 0, dates))
    xticks --> (xnums[inds], map(Dates.year, dates[inds]))

    # Shock contributions
    @series begin
        labels    = map(x -> x.name,  groups)
        cat_names = map(Symbol, labels)
        colors    = map(x -> x.color, groups)
        ngroups   = length(groups)

        labels     --> reshape(labels, 1, ngroups)
        color      --> reshape(colors, 1, ngroups)
        linealpha  --> 0
        bar_width  --> 1
        legendfont --> Plots.Font("sans-serif", 5, :hcenter, :vcenter, 0.0, colorant"black")

        inds = find(start_date .<= dates .<= end_date)
        x = df[inds, :date]
        y = convert(Array, df[inds, cat_names])
        StatPlots.GroupedBar((x, y))
    end

    seriestype := :line
    linewidth  := 2
    ylim       := :auto

    # Detrended mean history
    @series begin
        linecolor := hist_color
        label     := hist_label

        inds = intersect(find(start_date .<= dates .<= end_date),
                         find(hist.means[1, :date] .<= dates .<= hist.means[end, :date]))
        xnums[inds], df[inds, :detrendedMean]
    end

    # Detrended mean forecast
    @series begin
        linecolor := forecast_color
        label     := forecast_label

        inds = intersect(find(start_date .<= dates .<= end_date),
                         find(hist.means[end, :date] .<= dates .<= forecast.means[end, :date]))
        xnums[inds], df[inds, :detrendedMean]
    end
end


text:
"""
    MyWorkflow
This is an example of Julia module. These docstrings can be seen on github
"""
module MyWorkflow

export func

"""
    func(x)
Returns double the numer `x` plus `1`.

```jldoctest
julia> func(2)
5
```
"""
func(x) = 2x + 1

"""
    hello(msg::String)
print string `Hello` followed by `msg`

```jldoctest
julia> MyWorkflow.hello("world")
"helloworld"
```
"""
hello(msg::String) = "hello" * msg

end


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# camera/analog.jl: Camera object analog features

"""
  gain(::Camera) -> (Float,String)
  Camera exposure gain number (dB) and automatic gain mode.
"""
gain(cam::Camera) = (get(SpinFloatNode(cam, "Gain")), get(SpinEnumNode(cam, "GainAuto")))

"""
  gain!(::Camera)

  Set automatic exposure gain on camera.
"""
gain!(cam::Camera) = set!(SpinEnumNode(cam, "GainAuto"), "Continuous")

"""
  gain!(::Camera, ::Number) -> Float

  Set exposure gain on camera to specified number (dB). Gain is clamped to
  range supported by camera. This function disables automatic gain.
"""
function gain!(cam::Camera, g)
  set!(SpinEnumNode(cam, "GainAuto"), "Off")
  set!(SpinFloatNode(cam, "Gain"), g)
end

"""
  gain_limits(::Camera) -> (Float, Float)

  Camera gain limits in dB.
"""
gain_limits(cam::Camera) = range(SpinFloatNode(cam,"Gain"))

"""
  autogain_limits!(::Camera, (lower, upper); clampwarn=true) -> (Float, Float)

  Write lower and upper limits of the Auto Gain Time (us) value.
  Values exceeding range are clamped to the allowable extrema and returned, and
  a warning is given, which can be disabled with `clampwarn=false`.
"""
function autogain_limits!(cam::Camera, lims; clampwarn=true)
  set!(SpinFloatNode(cam, "AutoExposureGainLowerLimit"), lims[1], clampwarn=clampwarn)
  set!(SpinFloatNode(cam, "AutoExposureGainUpperLimit"), lims[2], clampwarn=clampwarn)
  autogain_limits(cam)
end

"""
  autogain_limits(::Camera) -> (Float, Float)

  Lower and upper limits of the Auto Gain Time (us) value.
"""
function autogain_limits(cam::Camera)
  (get(SpinFloatNode(cam, "AutoExposureGainLowerLimit")),
   get(SpinFloatNode(cam, "AutoExposureGainUpperLimit")))
end

"""
  gammaenable(::Camera) -> Bool

  Return status of gamma correction
"""
gammaenable(cam::Camera) = get(SpinBooleanNode(cam, "GammaEnable"))

"""
  gammaenable!(::Camera, ::Bool) -> Bool

  Enable or disable gamma correction on camera.
"""
gammaenable!(cam::Camera, en::Bool) = set!(SpinBooleanNode(cam, "GammaEnable"), en)


text:
function get_favorites_list(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/favorites/list.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
module Model

using Distributions
using StatsFuns
import Printf

# For plotting.
using StatsPlots, LaTeXStrings
import StatsPlots.KernelDensity.kde
import LinearAlgebra

include("../Util/Util.jl")
using .Util: SkewT

include("../MCMC/MCMC.jl")
include("Simulate.jl")
include("Data.jl")
include("Prior.jl")
include("State.jl")
include("Tuners.jl")
include("Flag.jl")
include("loglike.jl")
include("updates/update.jl")
include("fit.jl")

include("PriorAM.jl")
include("StateAM.jl")
include("logprioram.jl")
include("loglikeam.jl")
include("amfit.jl")

include("logprior.jl")
include("pseudo_prior_updates/update.jl")
include("fit_via_pseudo_prior.jl")
include("postprocess.jl")

export UpdateBetaWithSkewT, UpdateLambdaWithSkewT, fit, amfit
export fit_via_pseudo_prior, ppfit, posterior_prob1, bayes_factor, dic

end # module

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation