text:
using LinearAlgebra

############################
## TYPE DEFINITION ####
############################
abstract type PhylogeneticDistance end
mutable struct Hamming <: PhylogeneticDistance end
mutable struct Correlation <: PhylogeneticDistance end
mutable struct MutualInfo <: PhylogeneticDistance end
mutable struct pValue <: PhylogeneticDistance end
mutable struct mfDCA <: PhylogeneticDistance end
mutable struct plmDCA <: PhylogeneticDistance end
UnionDistances=Union{Hamming,Correlation,pValue,MutualInfo,mfDCA,plmDCA}

struct PhyloOut
	PhyloProfile::Array{Int8,2}
	list_domains::Array{String,1}
	list_species::Array{String,1}
	PhyloDistance::Array{Float64,2}
    result_unsort::Array{Any,2}
	result_sorted::Array{Any,2}
end

################################################
# METHODS DEFINITION
###############################################

function evaluate_distance(::Hamming,
                          P::Matrix{Int8})
	M,N=size(P)
	hamm_matrix=Matrix{Float64}(N,N)
	for i1=1:N
		for i2=(i1+1):N
			hamm_matrix[i1,i2]=sum(P[:,i1].!=P[:,i2])/N
		end
	end
	hamm_matrix+=hamm_matrix'
	return hamm_matrix
end

function evaluate_distance(::Correlation,
			   P::Matrix{Int8})

	return cor(P)
end

########################################
# MUTUAL INFORMATION
## to be implemented
#function evaluate_distance(::MutualInfo,
#                           P::Matrix{Int8})
#	P = ones(P)*2 - P
#    Z = Array{Int8}(P')
#	N, M = size(Z)
#	q = Int(maximum(Z))
#	theta = 0.2
#	pseudocount = 0.25 #for the 0,1 Ising model
#	#frequencies of one(s)
#	Pi_true, Pij_true, Meff, _ = compute_new_frequencies(Z, theta)
#	Pi, Pij = add_pseudocount(Pi_true, Pij_true, Meff, Float64(pseudocount))
#
#
#	#compute mutual information
#	mi = zeros(N,N)
#	for i = 1:N
#		for j = (i+1):N
#			Pij_00 = 1 + Pij[i,j] - Pi[i] - Pi[j]
#			Pij_01 = Pi[j] - Pij[i,j]
#			Pij_10 = Pi[i] - Pij[i,j]
#			Pij_11 = Pij[i,j]
#			Pi_0 = 1 - Pi[i]
#			Pj_0 = 1 - Pi[j]
#			Pi_1 = Pi[i]
#			Pj_1 = Pi[j]
#			mi[i,j] = Pij_00*log(Pij_00/Pi_0*Pj_0) + Pij_01*log(Pij_01/Pi_0*Pj_1)+Pij_10*log(Pij_10/Pi_1*Pj_0) + Pij_11*log(Pij_11/Pi_1*Pj_1)
#		end
#	end
#
#	mi += mi'
#	return mi
#end
################################

function evaluate_distance(::pValue,
                          P::Matrix{Int8})
	M,N = size(P)

	d_fisher = zeros(N,N)
	m11 = 0
	m12 = 0
	m21 = 0
	m22 = 0
	for i = 1:N
		for j = (i+1):N
			for lin = 1:M
				m11 = m11 + (Int(P[lin,i] == 0 && P[lin,j] == 0))
				m12 = m12 + (Int(P[lin,i] == 1 && P[lin,j] == 0))
				m21 = m21 + (Int(P[lin,i] == 0 && P[lin,j] == 1))
				m22 = m22 + (Int(P[lin,i] == 1 && P[lin,j] == 1))
			end
			res = FisherExactTest(m11,m12,m21,m22)
			pvalue_fisher = pvalue(res, tail=:both, method=:minlike)
			d_fisher[i,j] = pvalue_fisher
		end
	end
	d_fisher += d_fisher'
	return d_fisher

end

function evaluate_distance(::mfDCA,
                           P::Matrix{Int8};
                           pseudocount::Real=0.25, #for the 0,1 Ising model
                           theta::Float64=0.2)

	M, N = size(P)

	Pi_true, Pij_true, Meff, _ = compute_new_frequencies(P, theta)
	Pi, Pij = add_pseudocount(Pi_true, Pij_true, Meff, Float64(pseudocount))
	C = compute_C(Pi, Pij)

	mJ = -inv(cholesky(C))

	S = correct_APC(mJ)

    #remove diagonal elements (fields)
	S-=Diagonal(S)

	return S
end

##MY PLM DCA
#function evaluate_distance(::plmDCA,
#                           P::Matrix{Int8};
#                           lambdaJ::Real=0.02,
#                           lambdaH::Real=0.05,
#                           theta::Float64=0.2)
#
#    M,N=size(P)
#    W,Meff=compute_weights(P,theta)
#
#    #PlmAlg is a type
#	plmalg =PlmDCA.PlmAlg(method,verbose, epsconv ,maxit, boolmask)
#	plmvar = PlmDCA.PlmVar(N,M,q,q*q,gaugecol,lambdaJ,lambdaH,Z,W)
#    Jmat, pslike = PlmDCA.MinimizePLAsym(plmalg,plmvar)
#end
#

function evaluate_distance(::plmDCA,
                           P::Matrix{Int8})
#                          P::Array{Int8})

    lambdaJ::Real = 0.02
    lambdaH::Real = 0.05
    theta = 0.2
    epsconv::Real = 1.0e-5
    maxit::Int = 1000
    min_separation::Int = 1
    verbose::Bool = true
    method::Symbol = :LD_LBFGS
    boolmask::Union{Array{Bool,2},Nothing} = nothing
    gaugecol::Int = -1

    M,N = size(P)
#	P = ones(P)*2-P
	P = ones(size(P))*2-P

#   Z = Array{Int8}(P)
    Z = Array{Int8}(P')

	q = Int(maximum(Z))

	Pi_true, Pij_true, Meff, W = compute_new_frequencies(Z, theta)
#	plmalg =PlmDCA.PlmAlg(method,verbose, epsconv ,maxit, boolmask)
	plmalg =PlmDCA.PlmAlg(method,verbose, epsconv ,maxit)

#	plmvar = PlmDCA.PlmVar(N,M,q,q*q,gaugecol,lambdaJ,lambdaH,Z,W)
	plmvar = PlmDCA.PlmVar(N,M,q,q*q,lambdaJ,lambdaH,Z,W)

    Jmat, pslike = PlmDCA.MinimizePLAsym(plmalg,plmvar)

	J = switch_gauge(Jmat, plmvar, min_separation)

	S = correct_APC(J)
	#remove diagonal elements (fields...)
	S-= Diagonal(S)

	return S

end

function switch_gauge(Jmat::Array{Float64,2}, var, min_separation::Int)

    q = var.q
    N = var.N

    JJ = reshape(Jmat[1:end-q,:], q,q,N-1,N)
    Jtemp1 = zeros( q,q,Int(N*(N-1)/2))
    Jtemp2 = zeros( q,q,Int(N*(N-1)/2))

    l = 1

    for i = 1:(N-1)
        for j = (i+1):N
            Jtemp1[:,:,l] = JJ[:,:,j-1,i]; #J_ij as estimated from from g_i.
            Jtemp2[:,:,l] = JJ[:,:,i,j]'; #J_ij as estimated from from g_j.
            l=l+1;
        end
    end

    J1 = zeros(q,q,Int(N*(N-1)/2))
    J2 = zeros(q,q,Int(N*(N-1)/2))

    #switch to Ising gauge
    for l = 1:Int(N*(N-1)/2)
#        J1[:,:,l] = Jtemp1[:,:,l] - repmat(mean(Jtemp1[:,:,l],dims=1),q,1) -repmat(mean(Jtemp1[:,:,l],dims=2),1,q) .+ mean(Jtemp1[:,:,l])
#        J2[:,:,l] = Jtemp2[:,:,l] - repmat(mean(Jtemp2[:,:,l],dims=1),q,1) -repmat(mean(Jtemp2[:,:,l],dims=2),1,q) .+ mean(Jtemp2[:,:,l])
        J1[:,:,l] = Jtemp1[:,:,l] - repeat(mean(Jtemp1[:,:,l],dims=1),q,1) -repeat(mean(Jtemp1[:,:,l],dims=2),1,q) .+ mean(Jtemp1[:,:,l])
        J2[:,:,l] = Jtemp2[:,:,l] - repeat(mean(Jtemp2[:,:,l],dims=1),q,1) -repeat(mean(Jtemp2[:,:,l],dims=2),1,q) .+ mean(Jtemp2[:,:,l])
    end
    J = 0.5 * ( J1 + J2 )

    #switch to lattice-gas model
    lattice_gas_J = zeros(N,N)
    l = 1
    for i = 1:N-1
        for j = i+1:N
		lattice_gas_J[i,j] = J[1,1,l]-J[1,2,l]-J[2,1,l]+J[2,2,l]
		l += 1
        end
    end

    lattice_gas_J += lattice_gas_J'

    return lattice_gas_J
end

max_stars_repo_name: clayfos/PhyDCA

max_issues_repo_name: clayfos/PhyDCA

max_forks_repo_name: clayfos/PhyDCA


text:
"""
    struct SymSparseMatrixCSR{Bi,T,Ti<:Integer} <: AbstractSparseMatrix{T,Ti}
      uppertrian :: SparseMatrixCSR{Bi,T,Ti}
    end

Matrix type for storing symmetric sparse matrices in the
Compressed Sparse Row format with `Bi`-based indexing (typically 0 or 1).
Only the upper triangle is stored (including the non zero diagonal entries),
which is represented by a `SparseMatrixCSR`.
The standard way of constructing a `SymSparseMatrixCSR` is through the
[`symsparsecsr`](@ref) function.
"""
struct SymSparseMatrixCSR{Bi,T,Ti<:Integer} <: AbstractSparseMatrix{T,Ti}
    uppertrian :: SparseMatrixCSR{Bi,T,Ti}
end

"""
    symsparsecsr(args...;symmetrize::Bool=false)
    symsparsecsr(::Val{Bi},args...;symmetrize::Bool=false) where Bi

Create  a `SymSparseMatrixCSR` with `Bi`-based indexing (1 by default)
from the same `args...` as one constructs a `SparseMatrixCSC`
with the [`sparse`](@ref) function. If `symmetrize == false` (the default)
the given arguments should only describe the upper triangle
of the matrix (including non zero diagonal values). If `symmetrize == true`
a non symmetric input is accepted and it will be symmetrized in-place
(i.e., changing the input arguments).
"""
function symsparsecsr(::Val{Bi},I,J,V,args...;symmetrize::Bool=false) where Bi
  if symmetrize
    Tv = eltype(V)
    α = Tv(0.5)
    for k in 1:length(I)
      r = I[k]
      c = J[k]
      if r > c
        I[k] = c
        J[k] = r
      end
      if r != c
        V[k] = α*V[k]
      end
    end
  end
  SymSparseMatrixCSR(sparsecsr(Val(Bi),I,J,V,args...))
end
symsparsecsr(args...;kwargs...) = symsparsecsr(Val(1),args...;kwargs...)

size(A::SymSparseMatrixCSR) = size(A.uppertrian)
IndexStyle(::Type{<:SymSparseMatrixCSR}) = IndexCartesian()
function getindex(A::SymSparseMatrixCSR, x::Integer, y::Integer)
  getindex(A.uppertrian,min(x,y),max(x,y))
end
function setindex!(A::SymSparseMatrixCSR, v, x::Integer, y::Integer)
  setindex!(A.uppertrian,v,min(x,y),max(x,y))
end

getrowptr(S::SymSparseMatrixCSR) = getrowptr(S.uppertrian)
getnzval(S::SymSparseMatrixCSR) = getnzval(S.uppertrian)
getcolval(S::SymSparseMatrixCSR) = getcolval(S.uppertrian)

"""
    getBi(S::SymSparseMatrixCSR{Bi}) where {Bi}

Return `Bi`.
"""
getBi(S::SymSparseMatrixCSR{Bi}) where {Bi} = Bi

"""
    getoffset(S::SymSparseMatrixCSR{Bi}) where {Bi}

Return `1-Bi`. Useful to convert from 1-based to `Bi`-based indexing
(by subtracting the offset).
"""
getoffset(S::SymSparseMatrixCSR{Bi}) where Bi = getoffset(Bi)

"""
    issparse(S::SymSparseMatrixCSR)

Returns `true`.
"""
issparse(S::SymSparseMatrixCSR) = true

"""
    nnz(S::SymSparseMatrixCSR)

Returns the number of stored elements in a sparse array,
which correspond to the nonzero entries in the upper triangle and diagonal.
"""
nnz(S::SymSparseMatrixCSR) = nnz(S.uppertrian)

"""
    nonzeros(S::SymSparseMatrixCSR)

Return a vector (1-based) of the structural nonzero values in sparse array S.
This includes zeros that are explicitly stored in the sparse array,
which correspond to the nonzero entries in the upper triangle and diagonal.
The returned vector points directly to the internal nonzero storage of S,
and any modifications to the returned vector will mutate S as well.
"""
nonzeros(S::SymSparseMatrixCSR) = nonzeros(S.uppertrian)

"""
    colvals(S::SparseMatrixCSR)

Return a vector of the col indices of `S`. The stored values are indexes to arrays
with `Bi`-based indexing, but the `colvals(S)` array itself is a standard 1-based
Julia `Vector`.
Any modifications to the returned vector will mutate S as well.
Providing access to how the col indices are stored internally
can be useful in conjunction with iterating over structural
nonzero values. See also [`nonzeros`](@ref) and [`nzrange`](@ref).
"""
colvals(S::SymSparseMatrixCSR) = colvals(S.uppertrian)

"""
    nzrange(S::SymSparseMatrixCSR, row::Integer)

Return the range of indices to the structural nonzero values of a
sparse matrix row section being in the diagonal or upper triangle.
The returned range of indices is always 1-based even for `Bi != 1`.
"""
nzrange(S::SymSparseMatrixCSR, row::Integer) = nzrange(S.uppertrian, row)

"""
    findnz(S::SymSparseMatrixCSR)

Return a tuple `(I, J, V)` where `I` and `J` are the row and column 1-based indices
of the stored ("structurally non-zero in diagonal + upper trianle") values in sparse matrix A,
and V is a vector of the values. The returned vectors are newly allocated
and are unrelated to the internal storage of matrix `S`.
"""
findnz(S::SymSparseMatrixCSR) = findnz(S.uppertrian)

"""
    count(pred, S::SymSparseMatrixCSR)
    count(S::SymSparseMatrixCSR)

Count the number of elements in `nonzeros(S)` for which predicate `pred` returns `true`.
If  `pred` not given, it counts the number of `true` values.
"""
count(pred, S::SymSparseMatrixCSR) = count(pred, S.uppertrian)
count(S::SymSparseMatrixCSR) = count(i->true, S)

function LinearAlgebra.fillstored!(a::SymSparseMatrixCSR,v)
  LinearAlgebra.fillstored!(a.uppertrian,v)
  a
end

function LinearAlgebra.rmul!(a::SymSparseMatrixCSR,v::Number)
  LinearAlgebra.rmul!(a.uppertrian,v)
  a
end

function mul!(y::AbstractVector,A::SymSparseMatrixCSR,v::AbstractVector, α::Number, β::Number)
  A.uppertrian.n == size(v, 1) || throw(DimensionMismatch())
  A.uppertrian.m == size(y, 1) || throw(DimensionMismatch())
  if β != 1
    β != 0 ? rmul!(y, β) : fill!(y, zero(eltype(y)))
  end
  o = getoffset(A)
  for row = 1:size(y, 1)
    @inbounds for nz in nzrange(A,row)
      col = A.uppertrian.colval[nz]+o
      y[row] += A.uppertrian.nzval[nz]*v[col]*α
      row != col && (y[col] += A.uppertrian.nzval[nz]*v[row]*α)
    end
  end
  return y
end

function mul!(y::AbstractVector,A::SymSparseMatrixCSR,v::AbstractVector)
  A.uppertrian.n == size(v, 1) || throw(DimensionMismatch())
  A.uppertrian.m == size(y, 1) || throw(DimensionMismatch())

  fill!(y,zero(eltype(y)))
  o = getoffset(A)
  for row = 1:size(y, 1)
    @inbounds for nz in nzrange(A,row)
      col = A.uppertrian.colval[nz]+o
      y[row] += A.uppertrian.nzval[nz]*v[col]
      row != col && (y[col] += A.uppertrian.nzval[nz]*v[row])
    end
  end
  return y
end

*(A::SymSparseMatrixCSR, v::Vector) = (y = similar(v,size(A,1));mul!(y,A,v))

function show(io::IO, ::MIME"text/plain", S::SymSparseMatrixCSR)
    xnnz = nnz(S)
    print(io, S.uppertrian.m, "×", S.uppertrian.n, " ",
              typeof(S), " with ", xnnz, " stored ",
              xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        print(io, ":")
        show(IOContext(io, :typeinfo => eltype(S)), S)
    end
end
show(io::IO, S::SymSparseMatrixCSR) = show(io, S.uppertrian)

Base.convert(::Type{T},a::T) where T<:SymSparseMatrixCSR = a
function Base.convert(
  ::Type{SymSparseMatrixCSR{Bi,Tv,Ti}},a::SymSparseMatrixCSR) where {Bi,Tv,Ti}
  utrian = convert(SparseMatrixCSR{Bi,Tv,Ti},a.uppertrian)
  SymSparseMatrixCSR(utrian)
end

ext: jl


text:
using DynamicHMC, LogDensityProblems
using Distributions, Parameters, Random
using StatsPlots, PositiveFactorizations, Calculus
using Flux

include(string(pwd(),"/src/simulation.jl"))
include(string(pwd(),"/gali_bayesian.jl"))
include(string(pwd(),"/mcmc/priors_v2.jl"))
include(string(pwd(),"/mcmc/foos.jl"))

yy,shocks = simulate_dsge(GAMMA_0,GAMMA_1,PSI,PI,500)

prob = DSGE_Model(yy[:,2],1/3)

prob((bet = 0.99,epsilon = 6,theta=2/3,sig=1,s2=1,phi=1,phi_pi=1.5,phi_y=0.5/4,rho_v=0.5))

t = problem_transform(prob)

P = TransformedLogDensity(t,prob)

aa = [0.99,6,2/3,1,1,1,1.5,0.5/4,0.5]

true_vals = TransformVariables.inverse(t,(bet = 0.99,epsilon = 6,theta=2/3,sig=1,s2=1,phi=1,phi_pi=1.5,phi_y=0.5/4,rho_v=0.5))

LogDensityProblems.logdensity(P,true_vals)

#grads = Calculus.gradient(P)

function LogDensityProblems.capabilities(::TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:bet, :epsilon, :theta, :sig, :s2, :phi, :phi_pi, :phi_y, :rho_v),Tuple{TransformVariables.ScaledShiftedLogistic{Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}}},DSGE_Model{Array{Float64,1},Float64}})
    LogDensityProblems.LogDensityOrder{1}() # can do gradient
end

#LogDensityProblems.dimension(::NormalPosterior) = 2 # for this problem

function LogDensityProblems.logdensity_and_gradient(problem::TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:bet, :epsilon, :theta, :sig, :s2, :phi, :phi_pi, :phi_y, :rho_v),Tuple{TransformVariables.ScaledShiftedLogistic{Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}}},DSGE_Model{Array{Float64,1},Float64}}, x)
    logdens = LogDensityProblems.logdensity(P,x)
    grad = -1*Flux.gradient(y->model_flux(y)[1],x)[1] .*exp.(model_flux(x)[1])
    logdens, grad
end

res = @time mcmc_with_warmup(Random.GLOBAL_RNG, P,1000)

posterior = transform.(t, res.chain)

post_array = zeros(1000,9)

for i in 1:1000
    @unpack bet, epsilon, theta, sig, s2,phi,phi_pi,phi_y,rho_v = posterior[i,:][1]
    post_array[i,:] = [bet epsilon theta sig s2 phi phi_pi phi_y rho_v]
end

StatsPlots.density(post_array[:,1], legend = :none)
title!(latexify("beta()"))
vline!([true_pars[:bet]])
vline!([mean(post_array[:,1])])

StatsPlots.density(post_array[:,2], legend = :none)
title!(latexify("epsilon()"))
vline!([true_pars[:epsilon]])
vline!([mean(post_array[:,2])])

StatsPlots.density(post_array[:,3], legend = :none)
title!(latexify("theta()"))
vline!([true_pars[:theta]])
vline!([mean(post_array[:,3])])

StatsPlots.density(post_array[:,4], legend = :none)
title!(latexify("sigma()"))
vline!([true_pars[:sig]])
vline!([mean(post_array[:,4])])

StatsPlots.density(post_array[:,5], legend = :none)
title!(latexify("sigma()^2"))
vline!([true_pars[:s2]])
vline!([mean(post_array[:,5])])

StatsPlots.density(post_array[:,6], legend = :none)
title!(latexify("phi()"))
vline!([true_pars[:phi]])
vline!([mean(post_array[:,6])])

StatsPlots.density(post_array[:,7], legend = :none)
title!(latexify("phi()_pi",cdot = false))
vline!([true_pars[:phi_pi]])
vline!([mean(post_array[:,7])])

StatsPlots.density(post_array[:,8], legend = :none)
title!(latexify("phi()_y",cdot = false))
vline!([true_pars[:phi_y]])
vline!([mean(post_array[:,8])])

StatsPlots.density(post_array[:,9], legend = :none)
title!(latexify("rho()_v",cdot = false))
vline!([true_pars[:rho_v]])
vline!([mean(post_array[:,9])])

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
a=rand(Int,3).%9;b=rand(Int,3).%9;c=rand(Int,3).%9
abc=[Equation(:a,a),Equation(:b,b),Equation(:c,c)]

#Dot product
ex1=Ten(:a,:i)*Ten(:b,:i)
r1=ex1&abc
@assert r1==dot(a,b)

#Cross product
ex2=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)
r2=ex2&abc
@assert r2.x==cross(a,b)

ex3=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)*Ten(:b,:m)*Ten(:a,:m)+Ten(:c,:i)
r3=ex3&abc
@assert r3.x==cross(a,b)*dot(b,a)+c

#Both
ex4=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)*Ten(:a,:m)*Ten(:c,:m)
r4=ex4&abc
@assert r4.x==cross(a,b)*dot(a,c)

ex5=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)*Alt([:m,:n,:o])*Ten(:c,:n)*Ten(:b,:o)
r5=ex5&abc&Equation(:m,:i)
@assert r5==dot(cross(a,b),cross(c,b))

#Triangle area
ex6=Sqrt(Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)*Alt([:i,:m,:n])*Ten(:a,:m)*Ten(:b,:n))/2
r6=ex6&abc
@assert round(r6)==round(norm(cross(a,b))/2)

#Diadic products
ex7=Ten(:a,:i)*Ten(:b,:j)
r7=ex7&abc&@equs(i=2, j=3)
@assert r7==a[2]*b[3]

#Transpose
ex8=Equations.Transpose(Ten(a*b',[:i,:j]))
r8=simplify(ex8)
@assert r8.x==(a*b')'

#Determinant
A=rand(Int,3,3).%9
ex9a=Alt([:i,:j,:k])*Ten(A,[:i,1])*Ten(A,[:j,2])*Ten(A,[:k,3])
r9a=simplify(ex9a)
ex9b=Alt([:i,:j,:k])*Alt([:r,:s,:t])*Ten(A,[:i,:r])*Ten(A,[:j,:s])*Ten(A,[:k,:t])/6
r9b=simplify(ex9b)
@assert round(det(A))==r9a==round(r9b)

#Inverse
eq1=@equ detA=Alt([x,y,z])*Ten(A,[x,1])*Ten(A,[y,2])*Ten(A,[z,3])
eq2=@equ invA=Alt([j,m,n])*Alt([i,p,q])*Ten(A,[m,p])*Ten(A,[n,q])/(2detA)
eq3=eq2&eq1&Equation(:A,A)
@assert abs((eq3&@equs(i=1,j=3)).rhs-inv(A)[1,3])<0.001

#Matrix multiplication
B=rand(Int,3,3).%9
AB=[Equation(:A,A),Equation(:B,B)]
C=Ten(:A,[:i,:k])*Ten(:B,[:k,:j])
@assert C&AB&@equs(i=3,j=2) == (A*B)[3,2]

D=Ten(:A,[:i,:j])*Ten(:B,[:i,:j])
@assert D&AB==sum(A.*B)


text:
include("GB20.jl")
using MLBase
using DelimitedFiles
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV, Random
using DataFrames
using SharedArrays

rng = Random.seed!(2019)

gbdt = GB20.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                # Learning rate
  num_iterations = 100,               # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB20.GBLearner(gbdt, :regression)

############## σ² = 1 #################
Pval = []
for i in 1:1000
    if (10 <= i <= 19)
        push!(Pval, i)
    elseif (20 <= i <= 99)
        if (i % 10 == 0)
            push!(Pval, i)
        end
    else
        if (i % 100 == 0)
            push!(Pval, i)
        end

    end
end

Pval
temp_storedim = size(collect(Pval),1)

rmsefitted2 = []
HyperS = Array{Union{Missing, Nothing, Int, Float64}, 2}(nothing, temp_storedim, 1000)
TunedParams = Array{Union{Nothing, Float64}, 2}(nothing, temp_storedim, 3) # 28×3 Array, storing (α, τ, σ)

STORE = []
Pval
Pval
for pval in Pval
    # reproducibility
    rng = Random.seed!(2019)

    ### simulate the training and testing data
    X, Y = GB20.test_data(500; p = pval, sigma = 1.0)
    X_test, Y_test = GB20.test_data(500; p = pval, sigma = 1.0)
    μ = GB20.fried(X)
    μ_test = GB20.fried(X_test)

    hypers = GB20.Hypers(X, Y, 1, 1, 0.95, 0.5)
    hypers.σ = sqrt(1)
    hypers.τ = 3 / 10 * var(Y) / gbdt.num_iterations
    hypers.α = 0.95
    hypers.s = ones(size(X, 2)) / size(X, 2)

    # initialize opt
    opts = GB20.Opts()
    opts.num_burn = 5

    ### Exhaustive grid search for best param-combination ###
    param2update  = Dict(   :α => collect(range(0.8; length = 5, stop = 1.0)),
                            :τ => collect(range(0.05; length = 5, stop = 0.1)),
                            :σ => collect(range(1; length = 3, stop = 3)))
    BestHypers, BestCVError = GB20.TunePara_CV(X, Y, μ, gbdt, opts, hypers, 5, param2update; rng = rng)

    # store the tuned paramters:
    k = findall(x -> x == pval, Pval)
    TunedParams[k, :] = [BestHypers.α, BestHypers.τ, BestHypers.σ]
    #TunedParams[k, 1] = BestHypers.α
    #TunedParams[k, 2] = BestHypers.τ
    #TunedParams[k, 3] = BestHypers.σ

    ### Retrain the model with best params-combination and Evaluate on the Test set ###
    rmse2, Store_s = GB20.TBD_predict_and_rmse!(
        X,
        Y,
        X_test,
        Y_test,
        μ_test,
        gbdt,
        BestHypers,
        opts,
        rng)

    println("The rmse2 b.w. μ_test and pred for p = $(pval)  is $(rmse2).")

    push!(rmsefitted2, rmse2)
    push!(STORE, Store_s)

    # store the hyper.s for the ensembel tree
    j = findall(x -> x == pval, Pval)
    HyperS[j, 1:size(BestHypers.s, 1)] = BestHypers.s

end

rmsefitted2
HyperS
TunedParams

######## σ² = 10 #########
Pval = []
for i in 1:1000
    if (10 <= i <= 19)
        push!(Pval, i)
    elseif (20 <= i <= 99)
        if (i % 10 == 0)
            push!(Pval, i)
        end
    else
        if (i % 100 == 0)
            push!(Pval, i)
        end

    end
end

Pval
temp_storedim = size(collect(Pval),1)
rmsefitted2_sigma10 = []
HyperS_sigma10 = Array{Union{Missing, Nothing, Int, Float64}, 2}(nothing, temp_storedim, 1000)
TunedParams_sigma10 = Array{Union{Nothing, Float64}, 2}(nothing, temp_storedim, 3) # 28×3 Array, storing (α, τ, σ)

STORE_sigma10 = []
Pval

for pval in Pval[1]
    # reproducibility
    rng = Random.seed!(2019)

    ### simulate the training and testing data
    X, Y = GB20.test_data(500; p = pval, sigma = sqrt(10.0))
    X_test, Y_test = GB20.test_data(500; p = pval, sigma = sqrt(10.0))
    μ = GB20.fried(X)
    μ_test = GB20.fried(X_test)

    hypers = GB20.Hypers(X, Y, 1, 1, 0.95, 2.0)
    hypers.σ = sqrt(1)
    hypers.τ = 3 / 10 * var(Y) / gbdt.num_iterations
    hypers.α = 0.95
    hypers.s = ones(size(X, 2)) / size(X, 2)

    # initialize opt
    opts = GB20.Opts()
    opts.num_burn = 5

    ### Exhaustive grid search for best param-combination ###
    param2update  = Dict(   :α => collect(range(0.8; length = 5, stop = 1.0)),
                            :τ => collect(range(0.05; length = 5, stop = 0.1)),
                            :σ => collect(range(1; length = 3, stop = 3)))
    BestHypers, BestCVError = GB20.TunePara_CV(X, Y, μ, gbdt, opts, hypers, 5, param2update; rng = rng)

    # store the tuned paramters:
    k = findall(x -> x == pval, Pval)
    TunedParams_sigma10[k, :] = [BestHypers.α, BestHypers.τ, BestHypers.σ]
    #TunedParams[k, 1] = BestHypers.α
    #TunedParams[k, 2] = BestHypers.τ
    #TunedParams[k, 3] = BestHypers.σ

    # train the model and get the final model as well as rmse2
    rmse2, Store_s = GB20.TBD_predict_and_rmse!(
        X,
        Y,
        X_test,
        Y_test,
        μ_test,
        gbdt,
        BestHypers,
        opts,
        rng)

    println("The rmse2 b.w. μ_test and pred for p = $(pval)  is $(rmse2).")

    push!(rmsefitted2_sigma10, rmse2)
    push!(STORE_sigma10, Store_s)

    # store the hyper.s for the ensembel tree
    j = findall(x -> x == pval, Pval)
    HyperS_sigma10[j, 1:size(BestHypers.s, 1)] = BestHypers.s

    #println("After $(opts.num_burn) times of burning, the sum of the first 5 elements in the prob vec is: $(sum(hypers.s[1:5])).")
end

rmsefitted2_sigma10
HyperS_sigma10
TunedParams_sigma10

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
module CSwitch

export @cswitch

# inspired by dcjones's Switch.jl
macro cswitch(constexpr, body)
    case2label = Dict{Any,Symbol}()
    flow = Expr(:block)
    end_label = gensym("end")
    default_label = end_label

    for arg in body.args
        if Meta.isexpr(arg, :macrocall) && arg.args[1] == Symbol("@case")
            label = gensym("case")
            case2label[arg.args[3]] = label
            labelexpr = Expr(:symboliclabel, label)
            push!(flow.args, labelexpr)
        elseif Meta.isexpr(arg, :macrocall) && arg.args[1] == Symbol("@default")
            default_label = gensym("default")
            labelexpr = Expr(:symboliclabel, default_label)
            push!(flow.args, labelexpr)
        elseif arg == Expr(:break)
            labelexpr = Expr(:symbolicgoto, end_label)
            push!(flow.args, labelexpr)
        else
            push!(flow.args, arg)
        end
    end
    push!(flow.args, Expr(:symboliclabel, end_label))

    jumptable = Expr(:block)
    for (case, label) in case2label
        condition = Expr(:call, :(==), constexpr, case)
        push!(jumptable.args, Expr(:if, condition, Expr(:symbolicgoto, label)))
    end
    push!(jumptable.args[end].args, Expr(:symbolicgoto, default_label))

    return esc(Expr(:block, jumptable, flow))
end

end # module


text:
using ..timeevolution: nh_hamiltonian, dmaster_h!, dmaster_nh!, check_master

"""
    iterative!(rho0, H, J, [method!], args...; [log=false], kwargs...) -> rho[, log]

Compute the steady state density matrix of master equation defined by a
Hamiltonian and a set of jump operators by solving `L rho = 0` via an iterative
method provided as argument.

# Arguments
* `rho0`: Initial density matrix. Note that this gets mutated in-place.
* `H`: Operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators.
* `method!`: The iterative method to be used. Defaults to `IterativeSolvers.bicgstabl!`.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates) for
    the jump operators. If nothing is specified all rates are assumed to be 1.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
    operators. If they are not given they are calculated automatically.
* `args...`: Further arguments are passed on to the iterative solver.
* `kwargs...`: Further keyword arguments are passed on to the iterative solver.
See also: [`iterative`](@ref)

Credit for this implementation goes to Z. Denis and F. Vicentini.
See also https://github.com/Z-Denis/SteadyState.jl
"""
function iterative!(rho0::Operator, H::AbstractOperator, J,
                    method! = IterativeSolvers.bicgstabl!, args...;
                    rates=nothing, Jdagger=dagger.(J), kwargs...)

    # Solution x must satisfy L*x = y with y[end] = tr(x) = 1 and y[j≠end] = 0.
    M = length(rho0.basis_l)
    x0 = similar(rho0.data, M^2+1)
    x0[1:end-1] .= reshape(rho0.data, M^2)
    x0[end] = zero(eltype(rho0))

    y = similar(rho0.data, M^2+1)
    y[1:end-1] .= zero(eltype(rho0))
    y[end] = one(eltype(rho0))

    # Define the linear map lm: rho ↦ L(rho)
    lm = _linmap_liouvillian(rho0,H,J,Jdagger,rates)

    log = get(kwargs,:log,false)

    # Solve the linear system with the iterative solver, then devectorize rho
    if !log
        rho0.data .= @views reshape(method!(x0,lm,y,args...;kwargs...)[1:end-1],(M,M))
        return rho0
    else
        R, history = method!(x0,lm,y,args...;kwargs...)
        rho0.data .= @views reshape(R[1:end-1],(M,M))
        return rho0, history
    end
end

"""
    iterative(H, J, [method!], args...; [log=false], kwargs...) -> rho[, log]

Compute the steady state density matrix of master equation defined by a
Hamiltonian and a set of jump operators by solving `L rho = 0` via an iterative
method provided as argument.

# Arguments
* `rho0`: Initial density matrix. Note that this gets mutated in-place.
* `H`: Operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators.
* `method!`: The iterative method to be used. Defaults to `IterativeSolvers.bicgstabl!`.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates) for
    the jump operators. If nothing is specified all rates are assumed to be 1.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
    operators. If they are not given they are calculated automatically.
* `rho0=nothing`: Initial density operator.
* `args...`: Further arguments are passed on to the iterative solver.
* `kwargs...`: Further keyword arguments are passed on to the iterative solver.
See also: [`iterative!`](@ref)

Credit for this implementation goes to Z. Denis and F. Vicentini.
See also https://github.com/Z-Denis/SteadyState.jl
"""
function iterative(H::AbstractOperator, args...;
                rho0=nothing, kwargs...)
    if rho0 === nothing
        rho = DenseOperator(H.basis_l, H.basis_r)
        rho.data[1,1] = 1
    else
        rho = deepcopy(rho0)
    end
    return iterative!(rho, H, args...; kwargs...)
end

function _linmap_liouvillian(rho,H,J,Jdagger,rates)
    bl = rho.basis_l
    br = rho.basis_r
    M = length(bl)

    # Cache stuff
    drho = copy(rho)
    # rho = copy(rho)
    Jrho_cache = copy(rho)

    # Check reducibility
    isreducible = check_master(rho,H,J,Jdagger,rates)
    if isreducible
        Hnh = nh_hamiltonian(H,J,Jdagger,rates)
        Hnhdagger = dagger(Hnh)
        dmaster_ = (drho,rho) -> dmaster_nh!(drho,Hnh,Hnhdagger,J,Jdagger,rates,rho,Jrho_cache)
    else
        dmaster_ = (drho,rho) -> dmaster_h!(drho,H,J,Jdagger,rates,rho,Jrho_cache)
    end

    # Linear mapping
    function f!(y,x)
        # Reshape
        rho.data .= @views reshape(x[1:end-1], M, M)
        # Apply function
        dmaster_(drho,rho)
        # Recast data
        copyto!(y, 1, drho.data, 1, M^2)
        y[end] = tr(rho)
        return y
    end

    return LinearMaps.LinearMap{eltype(rho)}(f!,M^2+1;ismutating=true,issymmetric=false,ishermitian=false,isposdef=false)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions for debugging purposed.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    debug_keycode()

Debug key codes.

"""
function debug_keycode()
    # Initialize the terminal.
    term = REPL.Terminals.TTYTerminal("", stdin, stdout, stderr)

    # Switch the terminal to raw mode, meaning that all keystroke is immediatly
    # passed to us instead of waiting for <return>.
    REPL.Terminals.raw!(term, true)

    write(
        term.out_stream,
        "Type any key to echo the processed keycode. Hit q to exit.\n\n"
    )

    try
        while true
            k = _jlgetch(term.in_stream)
            println(k)
            k.value == "q" && break
        end
    finally
        REPL.Terminals.raw!(term, false)
    end

    return nothing
end


text:
struct IntervalLabels <: AbstractData end

Base.@kwdef mutable struct LabelledInterval <: AbstractModel
    label::String = ""
    nested_interval::NestedInterval
end

Base.@kwdef mutable struct LabelledIntervals{T₁ <: Dict{String, LabelledInterval}}  <: AbstractModel
    description::String = ""
    labelled_intervals::T₁
end

function get_description(lis::LabelledIntervals)
    lis.description
end

function get_labelled_intervals(lis::LabelledIntervals)
    lis.labelled_intervals
end

function get_label(li::LabelledInterval)
    li.label
end

function set_label!(li::LabelledInterval, label)
    li.label = label
end

function get_nested_interval(li::LabelledInterval)
    li.nested_interval
end

function set_nested_interval!(li::LabelledInterval, nested_interval)
    li.nested_interval = nested_interval
end

Base.@kwdef struct LabelIntervalDisplayProperties{T₁ <: NTuple} <: AbstractDisplayProperties
    id::String
    caption::String
    col::ImVec4 = ImVec4(0.9, 0, 0, 0.3)
    textcol::ImVec4 = ImVec4(0.4, 0.4, 0.4, 1)
    createwindow::Bool = true
    layout::RectangularLayout
    padding::T₁ = (0, 0, 0 ,0)
end

function get_id(p::LabelIntervalDisplayProperties)
    p.id
end

function get_caption(p::LabelIntervalDisplayProperties)
    p.caption
end

function get_col(p::LabelIntervalDisplayProperties)
    p.col
end

function get_textcol(p::LabelIntervalDisplayProperties)
    p.textcol
end

function is_new_window(p::LabelIntervalDisplayProperties)
    p.createwindow
end

function get_layout(p::LabelIntervalDisplayProperties)
    p.layout
end

function get_padding(p::LabelIntervalDisplayProperties)
    p.padding
end

function get_plotcontext(p::LabelIntervalDisplayProperties)
   p.plotcontext
end

mutable struct LabelIntervalControl <: AbstractControl
    isenabled::Bool
    nested_interval::NestedInterval
    buffer::String
end

struct LabelIntervalContext{T₁ <: AbstractControl,   T₂ <: AbstractModel,  T₃ <: AbstractDisplayProperties, T₄ <: AbstractPlotContext} <: AbstractContext
    control::T₁
    model::T₂
    display_properties::T₃
    plot_context::T₄
end

function (context::LabelIntervalContext{<: LabelIntervalControl,   <: LabelledIntervals,  <:LabelIntervalDisplayProperties, <: PlotContext})()
        control = context.control
        model = context.model
        display_properties = context.display_properties
        plot_context = context.plot_context
        id = get_id(display_properties)
        caption = get_caption(display_properties)
        is_new_window(display_properties) ? CImGui.Begin(caption*id,C_NULL, CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus) : nothing
        isenabled(control) ? control(model, display_properties, plot_context) : nothing
        is_new_window(display_properties) ? CImGui.End() : nothing
end

function (control::LabelIntervalControl)(model::LabelledIntervals, properties::LabelIntervalDisplayProperties, plot_context::PlotContext)
    id = get_id(properties)
    caption = get_caption(properties)
    col = get_col(properties)
    textcol = get_textcol(properties)
    rectangle = get_layout(properties)
    totalwidth = get_width(rectangle)
    totalheight = get_height(rectangle)
    padding = get_padding(properties)
    CImGui.Dummy(ImVec2(0, padding[1]))
    width = totalwidth - padding[2]
    height = totalheight
    CImGui.Indent(padding[2])
    pos = CImGui.GetCursorScreenPos()
    #CImGui.InvisibleButton("###Events Overview Button", ImVec2(width, height))
    CImGui.Unindent(padding[2])
    CImGui.SetCursorScreenPos(pos.x, pos.y)
    # Plot the overview
    plot_context()

    draw_list = CImGui.GetWindowDrawList()
    yoffset = 3
    # Draw the labelled intervals
    for (key, labelled_interval) in pairs(get_labelled_intervals(model))
        # Map interval position to screen coordinates
        nested_interval = get_nested_interval(labelled_interval)
        interval = get_interval(nested_interval)
        start = get_start(nested_interval)
        stop =  get_stop(nested_interval)
        x₀ = stretch_linearly(start, first(interval), last(interval), pos.x, pos.x + width)
        x₁ = stretch_linearly(stop, first(interval), last(interval), pos.x, pos.x + width)
        # Shaded region
        CImGui.AddRectFilled(draw_list, ImVec2(x₀, pos.y + yoffset), ImVec2(x₁,  pos.y + height + yoffset), Base.convert(ImU32, col))
        # Text Label
        first_nul = findfirst(isequal('\0'), key)
        str = isnothing(first_nul) ? key : key[1:first_nul-1]
        CImGui.AddText(draw_list, ImVec2(x₀ + 5, pos.y + yoffset + 5), Base.convert(ImU32, textcol) , str)
        # Handle event selection
        CImGui.SetCursorScreenPos(x₀, pos.y)
        CImGui.InvisibleButton("###interval button $str", ImVec2(x₁ - x₀, height))
        # TODO Improve tooltip
        #CImGui.IsItemHovered() && CImGui.SetTooltip(string("Hello"))
        if CImGui.IsItemClicked()
            # Update the principal nested interval in accordance with the selection.
            set_start!(control.nested_interval, start)
            set_stop!(control.nested_interval, stop)
            # Update the text for the label in accordance with the selection.
            buffer = key*"\0"^(1)
            # Allow up to 64 characters for the label buffer.
            pad_null = max(0, 255 - length(buffer) + 1)
            control.buffer = buffer*"\0"^(pad_null)
        end
        CImGui.SameLine()
        CImGui.SetCursorScreenPos(pos.x, pos.y)
    end

    # Draw indicators for the currently selected interval
    current_nested_interval = control.nested_interval
    interval = get_interval(current_nested_interval)
    start = get_start(current_nested_interval)
    stop =  get_stop(current_nested_interval)
    x₀ = stretch_linearly(start, first(interval), last(interval), pos.x, pos.x + width)
    x₁ = stretch_linearly(stop, first(interval), last(interval), pos.x, pos.x + width)

    CImGui.AddLine(draw_list, ImVec2(x₀, pos.y + yoffset), ImVec2(x₀ , pos.y + height + yoffset), Base.convert(ImU32, ImVec4(0.0, 0.0, 0.0, 0.9)), Cfloat(2));
    CImGui.AddLine(draw_list, ImVec2(x₁, pos.y + yoffset), ImVec2(x₁ , pos.y + height + yoffset), Base.convert(ImU32, ImVec4(0.0, 0.0, 0.0, 0.9)), Cfloat(2));

    # Handle creation of new event labels
    CImGui.NewLine()
    CImGui.NewLine()
    CImGui.Indent(padding[2])
    CImGui.PushItemWidth(150)
    CImGui.InputText("###event description", control.buffer, length(control.buffer))
    CImGui.PopItemWidth()
    CImGui.SameLine()
    CImGui.Button("Create Interval Label") && add_label!(get_labelled_intervals(model), control.buffer, control.nested_interval)
    CImGui.SameLine()
    CImGui.Button("Remove Interval Label") && remove_label!(get_labelled_intervals(model), control.buffer)
    CImGui.Unindent(padding[2])
end

function add_label!(labelled_intervals::Dict{String, LabelledInterval}, label::String, ni::NestedInterval)
    key = extract_string(label)
    labelled_intervals[key] = LabelledInterval(key, copy(ni))
end

function remove_label!(labelled_intervals::Dict{String, LabelledInterval}, label::String)
    key = extract_string(label)
    if haskey(labelled_intervals, key)
        delete!(labelled_intervals, key)
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashbio_oncoprint

"""
    dashbio_oncoprint(;kwargs...)

An OncoPrint component.
The OncoPrint component is used to view multiple genetic alteration events
through an interactive and zoomable heatmap. It is a React/Dash port of the
popular oncoPrint() function from the BioConductor R package.
Under the hood, the rendering is done using Plotly.js built upon D3.
Plotly's interactivity allows the user to bind clicks and hovers to genetic
events, allowing the user to create complex bioinformatic apps or workflows
that rely on crossfiltering.
Read more about the component here:
https://github.com/plotly/react-oncoprint
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique to the component.
- `backgroundcolor` (String; optional): Default color for the tracks, in common name, hex, rgb or rgba format.
If left blank, will default to a light grey rgb(190, 190, 190).
- `colorscale` (Bool | Dict; optional): If not null, will override the default OncoPrint colorscale.
Default OncoPrint colorscale same as CBioPortal implementation.
Make your own colrscale as a {'mutation': COLOR} dict.
Supported mutation keys are ['MISSENSE, 'INFRAME', 'FUSION',
'AMP', 'GAIN', 'HETLOSS', 'HMODEL', 'UP', 'DOWN']
Note that this is NOT a standard plotly colorscale.
- `data` (Array; optional): Input data, in CBioPortal format where each list entry is a dict
consisting of 'sample', 'gene', 'alteration', and 'type'
- `eventDatum` (Dict; optional): A Dash prop that returns data on clicking, hovering or resizing the viewer.
- `height` (Real | String; optional): Height of the OncoPrint.
Will disable auto-resizing of plots if set.
- `padding` (Real; optional): Adjusts the padding (as a proportion of whitespace) between two tracks.
Value is a ratio between 0 and 1.
Defaults to 0.05 (i.e., 5 percent). If set to 0, plot will look like a heatmap.
- `range` (Array; optional): .Toogles whether or not to show a legend on the right side of the plot,
with mutation information.
- `showlegend` (Bool; optional): .Toogles whether or not to show a legend on the right side of the plot,
with mutation information.
- `showoverview` (Bool; optional): .Toogles whether or not to show a heatmap overview of the tracks.
- `width` (Real | String; optional): Width of the OncoPrint.
Will disable auto-resizing of plots if set.
"""
function dashbio_oncoprint(; kwargs...)
        available_props = Symbol[:id, :backgroundcolor, :colorscale, :data, :eventDatum, :height, :padding, :range, :showlegend, :showoverview, :width]
        wild_props = Symbol[]
        return Component("dashbio_oncoprint", "OncoPrint", "dash_bio", available_props, wild_props; kwargs...)
end


text:
# This is an implementation of the CMSA-ES as described in the paper:
#  H. Beyer and B. Sendhoff, "Covariance Matrix Adaptation Revisited –
#  the CMSA Evolution Strategy", 2008
#

function normalize_utilities(utilities)
  utilities / sum(utilities)
end

function linear_utilities(mu, lambda)
  normalize_utilities(hcat(ones(1, mu), zeros(1, lambda-mu)))
end

function log_utilities(mu, lambda)
  normalize_utilities(hcat((log(mu+1) - log(1:mu))', zeros(1, lambda-mu)))
end

# Optimize the n-dimensional objective func with a (mu,lambda) CMSA-ES.
function cmsa_es(p;
  trace = true,

  # Stopping criteria related
  max_seconds = 4*numdims(p), max_evals_per_dim = 1e7,
  ftol = 1e-7, xtol = 1e-10, stol = 1e-10,
  max_rounds_without_improvement = 500,

  # Starting points, will be random unless specified
  xmean = false,

  # Algorithm specific params:
  covarMatrixSampler = CholeskyCovarSampler,
  utilitiesFunc = log_utilities,
  lambda = 4*numdims(p),
  mu = int(max(ceil(lambda/rand(4:20)), 1.0)),
  tau = 1 / sqrt(2*numdims(p)), # Equation (1) on page 5 in Beyer2008
  tau_c = 1 + numdims(p) * (numdims(p) + 1) / (2 * mu),  # Equation (2) on page 5 in Beyer2008
  sigma = 0.05*rand(1:8)*minimum(diameters(search_space(p))),
  decompose_covar_prob = 0.4

  )

  N = numdims(p)
  ss = search_space(p)
  max_evals = max_evals_per_dim * N
  a = 1 - 1 / tau_c
  C = covarMatrixSampler(N)
  utilities = utilitiesFunc(mu, lambda)
  xbest = xmean = rand_individual(ss)   # Current best mean value.
  fbest = eval1(xbest, p)
  num_fevals = 1
  archive = BlackBoxOptim.TopListArchive(N, 10)
  add_candidate!(archive, fbest, xbest[:], num_fevals)
  fevals_last_best = num_fevals
  next_print_covar = 100
  termination_reason = "?"

  start_time = time()

  # Now lets optimize! Ensure we run at least one iteration.
  while(true)

    if (time() - start_time) > max_seconds
      termination_reason = "Exceeded time budget"
      break
    end

    if num_fevals > max_evals
      termination_reason = "Exceeded function eval budget"
      break
    end

    # Decompose only with a given probability => saves time
    if rand() <= decompose_covar_prob
      decompose!(C)
    end

    # Generate new population
    sigmas = sigma * exp( tau * randn(1, lambda) )  # 1*lambda
    s = multivariate_normal_sample(C, N, lambda)
    z = broadcast(*, sigmas, s)                     # n*lambda
    xs = repmat(xmean, 1, lambda) + z               # n*lambda

    # Evaluate fitness
    fitnesses = eval_fitnesses(p, xs, lambda)
    num_fevals += lambda

    # Check if best new fitness is best ever and print some info if tracing.
    indbest = argmin(fitnesses)
    fbest_new = fitnesses[indbest]

    if fbest_new < fbest
      xbest = xs[:, indbest]
      fbest = fbest_new
      add_candidate!(archive, fbest, xbest, num_fevals)
      fevals_last_best = num_fevals

      if fitness_is_within_ftol(p, ftol, fbest)
        termination_reason = "Within ftol"
        break
      end

      if trace
        println("$(num_fevals): Best fitness = $(fbest)")
        if num_fevals > next_print_covar
          next_print_covar = num_fevals + 100
          println("covar summary: ", sumstats(C.C, (x) -> @sprintf("%.2e", x)))
          println("sigma: ", sigma)
        end
      end
    else
      if (num_fevals - fevals_last_best) > max_rounds_without_improvement * lambda
        termination_reason = "Max rounds without improvement reached"
        break
      end
    end

    # Assign weights to the best individuals according to the utilities vector.
    weights = assign_weights(lambda, fitnesses, utilities)

    # Calculate new mean value based on weighted version of steps in population.
    xmean += (z * weights)

    # Update the covariance matrix
    uc = zeros(Float64, N, N)
    for i in 1:lambda
      if weights[i] > 0.0
        se = s[:,i]
        uc += weights[i] * (se * se')
      end
    end
    update_covariance_matrix!(C, uc, a)
    #ws = broadcast(*, weights', s)
    #update_covariance_matrix!(C, (ws * ws'), covar_learning_rate)

    # Adapt sigma for next round
    sigma = sigmas * weights

    # Terminate if sigma very small
    #println("sigma = $(sigma[1,1])")
    if sigma[1,1] < stol
      termination_reason = "Sigma too small"
      break
    end
  end

  return xbest, fbest, num_fevals, termination_reason, archive
end

function eval_fitnesses(problem, xs, lambda = size(xs, 2))
  fitnesses = zeros(lambda)
  for i in 1:lambda
    fitnesses[i] = eval1(xs[:,i], problem)
  end
  fitnesses
end

function assign_weights(n, fitnesses, utilities; minimize = true)
  us_ordered = zeros(n, 1)
  perms = sortperm(fitnesses, rev = !minimize)
  for i in 1:n
    us_ordered[perms[i]] = utilities[i]
  end
  us_ordered
end

# We create different types of Covariance matrix samplers based on different
# decompositions.
abstract type CovarianceMatrixSampler end

function update_covariance_matrix!(cms::CovarianceMatrixSampler, delta, a)
  C = a * cms.C + (1 - a) * delta
  cms.C = triu(C) + triu(C,1)' # Ensure C is symmetric. Should not be needed, investigate...
end

mutable struct EigenCovarSampler <: CovarianceMatrixSampler
  C::Array{Float64,2}
  B::Array{Float64,2}
  diagD::Array{Float64,1}

  EigenCovarSampler(n) = begin
    new(Matrix{Float64}(I, n,n), Matrix{Float64}(I, n,n), ones(n))
  end
end

function decompose!(cms::EigenCovarSampler)
  try
    EV, B = eig(cms.C)
    cms.B = B
    cms.diagD = sqrt(EV)
  catch
    # We don't update if there is some problem
  end
end

function multivariate_normal_sample(cms::CovarianceMatrixSampler, n, m)
  cms.B * (cms.diagD .* randn(n, m))
end

mutable struct CholeskyCovarSampler <: CovarianceMatrixSampler
  C::Array{Float64,2}
  sqrtC::Array{Float64,2}

  CholeskyCovarSampler(n) = begin
    new(Matrix{Float64}(I, n,n), Matrix{Float64}(I, n,n))
  end
end

function decompose!(cms::CholeskyCovarSampler)
  try
    cms.sqrtC = chol(cms.C)'
  catch error
    # We don't update if there is some problem
  end
end

function multivariate_normal_sample(cms::CholeskyCovarSampler, n, m)
  cms.sqrtC * randn(n, m)
end


text:
@testset "Maximum weighted matching" begin
    function neighbors(sym_matrix, i)
        return findall(v -> v != 0, sym_matrix[:, i])
    end

    @testset "Bipartite example" begin

        # from slide 65 of https://www.slideshare.net/KuoE0/acmicpc-bipartite-matching
        # but set 1 <-> 4 to 4.9 to avoid matching 1 <-> 4 and 3 <-> 5 which also sums to 10
        weight_matrix = [
            0 0 0 4.9 0 2
            0 0 0 3 1 -2
            0 0 0 0 5 0
            4.9 3 5 0 0 0
            0 1 0 0 0 0
            2 -2 0 0 0 0
        ]
        n = 6

        m = Model(CSJuMPTestOptimizer())
        @variable(m, x[1:n, 1:n], Bin)
        for i in 1:n, j in 1:n
            if weight_matrix[i, j] == 0 || i > j
                @constraint(m, x[i, j] == 0)
            end
        end
        for i in 1:n
            @constraint(
                m,
                sum(x[min(i, j), max(i, j)] for j in neighbors(weight_matrix, i)) <= 1
            )
        end

        @objective(m, Max, sum(weight_matrix .* x))
        optimize!(m)
        @test JuMP.objective_value(m) ≈ 10
        @test JuMP.value(x[1, 6]) == 1
        @test JuMP.value(x[2, 4]) == 1
        @test JuMP.value(x[3, 5]) == 1
        @test sum(JuMP.value.(x[1, 1:6])) == 1
        @test sum(JuMP.value.(x[2, 1:6])) == 1
        @test sum(JuMP.value.(x[3, 1:6])) == 1
        @test sum(JuMP.value.(x[4, 1:6])) == 0
        @test sum(JuMP.value.(x[5, 1:6])) == 0
        @test sum(JuMP.value.(x[6, 1:6])) == 0
    end
end


text:
module v2_tiny
#Tiny Yolo V2-tiny model configuration

using Knet: Knet, progress, progress!, gpu, KnetArray, relu, minibatch, conv4, pool, softmax

import ..sigmoid, ..Settings

mutable struct v2TinyChain
    layers
    v2TinyChain(layers...) = new(layers)
end

mutable struct Conv #Define convolutional layer
    w
    b
    stride
    padding
    f
end
mutable struct YoloPad #Define Yolo padding layer (assymetric padding).
    w
    grid_x::Int
    grid_y::Int
    minibatch_size::Int
end
struct Pool # Define pool layer
    size
    stride
    pad
end

YoloPad(w1::Int, w2::Int, cx::Int, cy::Int, grid_x::Int, grid_y::Int, minibatch_size::Int) =
    YoloPad(zeros(Float32, w1, w2, cx, cy), grid_x, grid_y, minibatch_size)#Constructor for Yolopad

Conv(w1::Int, w2::Int, cx::Int, cy::Int, st, pd, f) =
    Conv(randn(Float32, w1, w2, cx, cy), randn(Float32, 1, 1, cy, 1), st, pd, f)#Constructor for convolutional layer

#Assymetric padding function
function (y::YoloPad)(x)
    x = reshape(x, y.grid_x + 1, y.grid_y + 1, 1, 512 * y.minibatch_size)
    return reshape(conv4(y.w, x; stride = 1), y.grid_x, y.grid_y, 512, y.minibatch_size)
end

(p::Pool)(x) = pool(x; window = p.size, stride = p.stride, padding = p.pad) #pool function

(c::v2TinyChain)(x) = (for l in c.layers #chain function
    x = l(x)
end; x)

(c::Conv)(x) =
    c.f.(conv4(c.w, x; stride = c.stride, padding = c.padding) .+ c.b) #convolutional layer function

#leaky function
function leaky(x)
    return max(Float32(0.1) * x, x)
end

"""
    load(sets::Settings)

Loads the YOLOv2_tiny model.

Note on Asymmetric Padding (see https://github.com/Ybakman/YoloV2#asymmetric-padding)
An alternative solution for asymmetric padding on pooling. We solved the asymmetric padding process in the following 4 steps: Let's assume we have a matrix as (d1,d2,depth,minibatch_size)

1- Apply symmetric padding = 1 with the pooling function

2- reshape it into => (d1,d2,1,minibatch_size*depth)

3- apply the following convolutional layer => 1 0  with stride = 1, padding = 0
                                              0 0

4- Lastly, reshape again into => (d1,d2,depth,minibatch_size)

Now, we made pooling with asymmetric padding on the right and bottom side. If you want to
make it on left and top side, apply this conv layer =>  0 0
                                                        0 1

The last 3 steps are implemented in YoloPad(x) function
The symmetric padding part is implemented in the pooling function just before Yolopad(x)
"""
function load(sets::Settings)
return v2TinyChain(
    Conv(3, 3, 3, 16, 1, 1, leaky),
    Pool(2, 2, 0),
    Conv(3, 3, 16, 32, 1, 1, leaky),
    Pool(2, 2, 0),
    Conv(3, 3, 32, 64, 1, 1, leaky),
    Pool(2, 2, 0),
    Conv(3, 3, 64, 128, 1, 1, leaky),
    Pool(2, 2, 0),
    Conv(3, 3, 128, 256, 1, 1, leaky),
    Pool(2, 2, 0),
    Conv(3, 3, 256, 512, 1, 1, leaky),
    Pool(2, 1, 1),
    YoloPad(2, 2, 1, 1, sets.grid_x, sets.grid_y, sets.minibatch_size),
    Conv(3, 3, 512, 1024, 1, 1, leaky),
    Conv(3, 3, 1024, 1024, 1, 1, leaky),
    Conv(1, 1, 1024, sets.cell_bboxes * (5 + sets.num_classes), 1, 0, identity),
)
end

include(joinpath(@__DIR__, "loadweights.jl"))

end #module


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

module Collections

import Base: setindex!, done, get, hash, haskey, isempty, length, next, getindex, start, copymutable
import ..Order: Forward, Ordering, lt

export
    PriorityQueue,
    dequeue!,
    enqueue!,
    heapify!,
    heapify,
    heappop!,
    heappush!,
    isheap,
    peek

# Some algorithms that can be defined only after infrastructure is in place
Base.append!(a::Vector, iter) = _append!(a, Base.iteratorsize(iter), iter)

function _append!(a, ::Base.HasLength, iter)
    n = length(a)
    resize!(a, n+length(iter))
    @inbounds for (i,item) in zip(n+1:length(a), iter)
        a[i] = item
    end
    a
end

function _append!(a, ::Base.IteratorSize, iter)
    for item in iter
        push!(a, item)
    end
    a
end

# Heap operations on flat arrays
# ------------------------------

# Binary heap indexing
heapleft(i::Integer) = 2i
heapright(i::Integer) = 2i + 1
heapparent(i::Integer) = div(i, 2)

# Binary min-heap percolate down.
function percolate_down!(xs::AbstractArray, i::Integer, x=xs[i], o::Ordering=Forward, len::Integer=length(xs))
    @inbounds while (l = heapleft(i)) <= len
        r = heapright(i)
        j = r > len || lt(o, xs[l], xs[r]) ? l : r
        if lt(o, xs[j], x)
            xs[i] = xs[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
end

percolate_down!(xs::AbstractArray, i::Integer, o::Ordering, len::Integer=length(xs)) = percolate_down!(xs, i, xs[i], o, len)

# Binary min-heap percolate up.
function percolate_up!(xs::AbstractArray, i::Integer, x=xs[i], o::Ordering=Forward)
    @inbounds while (j = heapparent(i)) >= 1
        if lt(o, x, xs[j])
            xs[i] = xs[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
end

percolate_up!{T}(xs::AbstractArray{T}, i::Integer, o::Ordering) = percolate_up!(xs, i, xs[i], o)

"""
    heappop!(v, [ord])

Given a binary heap-ordered array, remove and return the lowest ordered element.
For efficiency, this function does not check that the array is indeed heap-ordered.
"""
function heappop!(xs::AbstractArray, o::Ordering=Forward)
    x = xs[1]
    y = pop!(xs)
    if !isempty(xs)
        percolate_down!(xs, 1, y, o)
    end
    x
end

"""
    heappush!(v, x, [ord])

Given a binary heap-ordered array, push a new element `x`, preserving the heap property.
For efficiency, this function does not check that the array is indeed heap-ordered.
"""
function heappush!(xs::AbstractArray, x, o::Ordering=Forward)
    push!(xs, x)
    percolate_up!(xs, length(xs), x, o)
    xs
end

# Turn an arbitrary array into a binary min-heap in linear time.
"""
    heapify!(v, ord::Ordering=Forward)

In-place [`heapify`](:func:`heapify`).
"""
function heapify!(xs::AbstractArray, o::Ordering=Forward)
    for i in heapparent(length(xs)):-1:1
        percolate_down!(xs, i, o)
    end
    xs
end

"""
    heapify(v, ord::Ordering=Forward)

Returns a new vector in binary heap order, optionally using the given ordering.
```jldoctest
julia> a = [1,3,4,5,2];

julia> Base.Collections.heapify(a)
5-element Array{Int64,1}:
 1
 2
 4
 5
 3

julia> Base.Collections.heapify(a, Base.Order.Reverse)
5-element Array{Int64,1}:
 5
 3
 4
 1
 2
```
"""
heapify(xs::AbstractArray, o::Ordering=Forward) = heapify!(copymutable(xs), o)

"""
    isheap(v, ord::Ordering=Forward)

Return `true` if an array is heap-ordered according to the given order.

```jldoctest
julia> a = [1,2,3]
3-element Array{Int64,1}:
 1
 2
 3

julia> Base.Collections.isheap(a,Base.Order.Forward)
true

julia> Base.Collections.isheap(a,Base.Order.Reverse)
false
```
"""
function isheap(xs::AbstractArray, o::Ordering=Forward)
    for i in 1:div(length(xs), 2)
        if lt(o, xs[heapleft(i)], xs[i]) ||
           (heapright(i) <= length(xs) && lt(o, xs[heapright(i)], xs[i]))
            return false
        end
    end
    true
end

# PriorityQueue
# -------------

"""
    PriorityQueue(K, V, [ord])

Construct a new [`PriorityQueue`](:obj:`PriorityQueue`), with keys of type
`K` and values/priorites of type `V`.
If an order is not given, the priority queue is min-ordered using
the default comparison for `V`.

A `PriorityQueue` acts like a `Dict`, mapping values to their
priorities, with the addition of a `dequeue!` function to remove the
lowest priority element.

```jldoctest
julia> a = Base.Collections.PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
Base.Collections.PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2
```
"""
type PriorityQueue{K,V,O<:Ordering} <: Associative{K,V}
    # Binary heap of (element, priority) pairs.
    xs::Array{Pair{K,V}, 1}
    o::O

    # Map elements to their index in xs
    index::Dict{K, Int}

    function PriorityQueue(o::O)
        new(Array{Pair{K,V}}(0), o, Dict{K, Int}())
    end

    PriorityQueue() = PriorityQueue{K,V,O}(Forward)

    function PriorityQueue(ks::AbstractArray{K}, vs::AbstractArray{V},
                           o::O)
        # TODO: maybe deprecate
        if length(ks) != length(vs)
            throw(ArgumentError("key and value arrays must have equal lengths"))
        end
        PriorityQueue{K,V,O}(zip(ks, vs), o)
    end

    function PriorityQueue(itr, o::O)
        xs = Array{Pair{K,V}}(length(itr))
        index = Dict{K, Int}()
        for (i, (k, v)) in enumerate(itr)
            xs[i] = Pair{K,V}(k, v)
            if haskey(index, k)
                throw(ArgumentError("PriorityQueue keys must be unique"))
            end
            index[k] = i
        end
        pq = new(xs, o, index)

        # heapify
        for i in heapparent(length(pq.xs)):-1:1
            percolate_down!(pq, i)
        end

        pq
    end
end

PriorityQueue(o::Ordering=Forward) = PriorityQueue{Any,Any,typeof(o)}(o)
PriorityQueue{K,V}(::Type{K}, ::Type{V}, o::Ordering=Forward) = PriorityQueue{K,V,typeof(o)}(o)

# TODO: maybe deprecate
PriorityQueue{K,V}(ks::AbstractArray{K}, vs::AbstractArray{V},
                   o::Ordering=Forward) = PriorityQueue{K,V,typeof(o)}(ks, vs, o)

PriorityQueue{K,V}(kvs::Associative{K,V}, o::Ordering=Forward) = PriorityQueue{K,V,typeof(o)}(kvs, o)

PriorityQueue{K,V}(a::AbstractArray{Tuple{K,V}}, o::Ordering=Forward) = PriorityQueue{K,V,typeof(o)}(a, o)

length(pq::PriorityQueue) = length(pq.xs)
isempty(pq::PriorityQueue) = isempty(pq.xs)
haskey(pq::PriorityQueue, key) = haskey(pq.index, key)

"""
    peek(pq)

Return the lowest priority key from a priority queue without removing that
key from the queue.
"""
peek(pq::PriorityQueue) = pq.xs[1]

function percolate_down!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while (l = heapleft(i)) <= length(pq)
        r = heapright(i)
        j = r > length(pq) || lt(pq.o, pq.xs[l].second, pq.xs[r].second) ? l : r
        if lt(pq.o, pq.xs[j].second, x.second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

function percolate_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        if lt(pq.o, x.second, pq.xs[j].second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

# Equivalent to percolate_up! with an element having lower priority than any other
function force_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        pq.index[pq.xs[j].first] = i
        pq.xs[i] = pq.xs[j]
        i = j
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

function getindex{K,V}(pq::PriorityQueue{K,V}, key)
    pq.xs[pq.index[key]].second
end

function get{K,V}(pq::PriorityQueue{K,V}, key, deflt)
    i = get(pq.index, key, 0)
    i == 0 ? deflt : pq.xs[i].second
end

# Change the priority of an existing element, or equeue it if it isn't present.
function setindex!{K,V}(pq::PriorityQueue{K, V}, value, key)
    if haskey(pq, key)
        i = pq.index[key]
        oldvalue = pq.xs[i].second
        pq.xs[i] = Pair{K,V}(key, value)
        if lt(pq.o, oldvalue, value)
            percolate_down!(pq, i)
        else
            percolate_up!(pq, i)
        end
    else
        enqueue!(pq, key, value)
    end
    value
end

"""
    enqueue!(pq, k, v)

Insert the a key `k` into a priority queue `pq` with priority `v`.

```jldoctest
julia> a = Base.Collections.PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
Base.Collections.PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2

julia> Base.Collections.enqueue!(a, "d", 4)
Base.Collections.PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 4 entries:
  "c" => 1
  "b" => 3
  "a" => 2
  "d" => 4
```
"""
function enqueue!{K,V}(pq::PriorityQueue{K,V}, key, value)
    if haskey(pq, key)
        throw(ArgumentError("PriorityQueue keys must be unique"))
    end
    push!(pq.xs, Pair{K,V}(key, value))
    pq.index[key] = length(pq)
    percolate_up!(pq, length(pq))
    pq
end

"""
    dequeue!(pq)

Remove and return the lowest priority key from a priority queue.

```jldoctest
julia> a = Base.Collections.PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
Base.Collections.PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2

julia> Base.Collections.dequeue!(a)
"c"

julia> a
Base.Collections.PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 2 entries:
  "b" => 3
  "a" => 2
```
"""
function dequeue!(pq::PriorityQueue)
    x = pq.xs[1]
    y = pop!(pq.xs)
    if !isempty(pq)
        pq.xs[1] = y
        pq.index[y.first] = 1
        percolate_down!(pq, 1)
    end
    delete!(pq.index, x.first)
    x.first
end

function dequeue!(pq::PriorityQueue, key)
    idx = pq.index[key]
    force_up!(pq, idx)
    dequeue!(pq)
    key
end

# Unordered iteration through key value pairs in a PriorityQueue
start(pq::PriorityQueue) = start(pq.index)

done(pq::PriorityQueue, i) = done(pq.index, i)

function next{K,V}(pq::PriorityQueue{K,V}, i)
    (k, idx), i = next(pq.index, i)
    return (pq.xs[idx], i)
end

end # module Collections

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using DifferentiableNAS
using Flux
using Flux: throttle, logitcrossentropy, onecold, onehotbatch, Optimiser
using Zygote: @nograd
using StatsBase: mean
using Parameters
using CUDA
using Distributions
using BSON
using Dates

include("../CIFAR10.jl")
include("../training_utils.jl")

argparams = trial_params()

num_ops = length(PRIMITIVES)

m = DARTSModelBN(num_cells = argparams.num_cells, channels = argparams.channels) |> gpu

function (hist::historiessml)()
    @show losses
    push!(hist.normal_αs_sm, m.normal_αs |> cpu)
    push!(hist.reduce_αs_sm, m.reduce_αs |> cpu)
    #push!(hist.activations, copy(m.activations.currentacts) |> cpu)
    push!(hist.train_losses, losses[1])
    push!(hist.val_losses, losses[2])
    CUDA.reclaim()
    GC.gc()
end

optimiser_α = Optimiser(WeightDecay(1e-3),ADAM(3e-4,(0.5,0.999)))
optimiser_w = Optimiser(WeightDecay(3e-4),Momentum(0.025, 0.9))

train, val = get_processed_data(argparams.val_split, argparams.batchsize, argparams.trainval_fraction)
test = get_test_data(argparams.test_fraction)

histepoch = historiessml()
histbatch = historiessml()
losses = [0.0, 0.0]

base_folder = prepare_folder("bndarts")

cbepoch = CbAll(CUDA.reclaim, histepoch, save_progress, CUDA.reclaim)
cbbatch = CbAll(CUDA.reclaim, histbatch, CUDA.reclaim)

for epoch in 1:argparams.epochs
    @show epoch
    DARTStrain1st!(loss, m, train, val, optimiser_α, optimiser_w, losses; cbepoch = cbepoch, cbbatch = cbbatch)
end

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


text:
function b2mc!(mc::AbstractVector, α)
    m = length(mc)
    o = zero(eltype(mc))

    mc[m] = mc[m]
    d = mc[m]
    for i=m-1:-1:1
        o = mc[i] + α*d
        d = mc[i]
        mc[i] = o
    end

    mc
end

b2mc(b::AbstractVector, α=0.35) = b2mc!(copy(b), α)

function b2mc!{T<:MelCepstrum}(state::SpectralParamState{T})
    assert_ready_to_filt(state)

    def = paramdef(state)
    α = allpass_alpha(def)
    data = rawdata(state)
    state.data = b2mc!(data, α)
    state.ready_to_filt = false
    state
end

function b2mc!{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T})
    throw(ArgumentError("""
                        filter coefficients based on mel-generalized cesptrum cannot be converted to mel-cepstrum"""))

end

function b2mc{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T})
    b2mc!(copy(state))
end


text:
import Plots: plot, xlabel!, ylabel!, title!
using DataFrames
using Gadfly

"Plot pY vs. ligand proportion."
function plotLigProp(rr, tps)
    GasProp = LinRange(0.0, 1.0, 10)
    pY = Array{Float64}(undef, length(tps), length(GasProp))

    for i = 1:length(GasProp)
        gas = GasProp[i]
        pros = 1.0 - gas
        totalconc = 70
        pYdata = TAMode.runTAM(tps, rr, (gas * totalconc, pros * totalconc))
        pY[:, i] = pYdata * TAMode.pYLS
    end

    pY1 = [[pY[1, :], pY[2, :]]]
    plot(GasProp, pY1, label = ["1 hr" "4 hr"], title = "Phosphorylated Receptor vs. Ligand Proportion", lw = 3)
    xlabel!("Proportion of Gas6")
    ylabel!("Phosphorylated Receptor (pY)")
end

"Plot Gas6 dose response."
function plotGasDose(rr, tps, conc; Gas = true)
    pY = Array{Float64}(undef, length(tps), length(conc))

    for i = 1:length(conc)
        if Gas
            stim = (conc[i], 0.0)
        else
            stim = (0.0, conc[i])
        end

        pY[:, i] = TAMode.runTAM(tps, rr, stim) * TAMode.pYLS
    end

    pY2g = [[pY[1, :], pY[2, :]]]
    plot(conc, pY2g, label = ["1 hr" "4 hr"], title = "Gas6 Dose Response", lw = 3)
    xlabel!("Gas6 Concentration")
    ylabel!("Phosphorylated Receptor (pY)")
end

"Plot dimer formation at 1 hr and 4 hr."
function plotDimers(rr, tps)
    GasProp = LinRange(0.0, 1.0, 10)
    GGdimers = vcat(zeros(9), ones(2), zeros(3), zeros(9), ones(2), zeros(5))
    PPdimers = vcat(zeros(11), ones(2), zeros(12), ones(2), zeros(3))
    GPdimers = vcat(zeros(13), 1, zeros(13), 1, zeros(2))
    GG = Array{Float64}(undef, length(tps), length(GasProp))
    PP = Array{Float64}(undef, length(tps), length(GasProp))
    GP = Array{Float64}(undef, length(tps), length(GasProp))

    for i = 1:length(GasProp)
        gas = GasProp[i]
        pros = 1.0 - gas
        totalconc = 70
        data = TAMode.runTAM(tps, rr, (gas * totalconc, pros * totalconc))
        GG[:, i] = data * GGdimers
        PP[:, i] = data * PPdimers
        GP[:, i] = data * GPdimers
    end

    plot1hr = [[GG[1, :], PP[1, :]], GP[1, :]]
    plot(GasProp, plot1hr, label = ["GG" "PP" "GP"], lw = 3)
    title!("Dimer Formation at 1 hr")
    xlabel!("Proportion of Gas6")
    ylabel!("Dimer Concentration")

    plot4hr = [[GG[2, :], PP[2, :]], GP[2, :]]
    plot(GasProp, plot4hr, label = ["GG" "PP" "GP"], lw = 3)
    title!("Dimer Formation at 4 hr")
    xlabel!("Proportion of Gas6")
    ylabel!("Dimer Concentration")
end

"Plot pY receptor expression"
function plotpYExpression(chn, AXLexpr, MerTKexpr, Tyro3expr, cellName)

    tps = range(0.0, 1440, length = 200)

    pYAXL = TAMode.pY .* TAMode.recpSpecific[1]
    pYTyro3 = TAMode.pY .* TAMode.recpSpecific[3]

    data = calcData(chn, AXLexpr, MerTKexpr, Tyro3expr, tps)
    pYA = (data * pYAXL)  # normal pYA
    pYT = (data * pYTyro3)  # normal pYT

    # AXL knockdown
    A_data = calcData(chn, 0.0, MerTKexpr, Tyro3expr, tps)
    A_pYA = (A_data * pYAXL) # AXLkd pYA
    A_pYT = (A_data * pYTyro3) # AXLkd pYT

    # Tyro3 knockdown
    T_data = calcData(chn, AXLexpr, MerTKexpr, 0.0, tps)
    T_pYA = (T_data * pYAXL) # Tyrokd pYA
    T_pYT = (T_data * pYTyro3) # Tyrokd pYT

    pyAXL = vcat(pYA[1], pYA[2], A_pYA[1], A_pYA[2], T_pYA[1], T_pYA[2])
    pyTyro3 = vcat(pYT[1], pYT[2], A_pYT[1], A_pYT[2], T_pYT[1], T_pYT[2])

    count = vcat(pyAXL, pyTyro3)
    treatment = repeat(["Normal", "Normal", "AXL Knockdown", "AXL Knockdown", "Tyro3 Knockdown", "Tyro3 Knockdown"], outer = 2)
    time = repeat(["1 hr", "8 hr"], outer = 6)
    receptorType = repeat(["Phosphyrlated AXL", "Phosphyrlated Tyro3"], inner = 6)

    D = DataFrame(Treatment = treatment, Time = time, ReceptorType = receptorType, Count = count)

    p2 = Gadfly.plot(
        D,
        x = :Time,
        y = :Count,
        color = :ReceptorType,
        xgroup = :Treatment,
        Geom.subplot_grid(Geom.bar(position = :stack)),
        Guide.xlabel("Treatment"),
        Guide.title(cellName),
    )

end

"Plot expression over time"
function plotTimeSeries(chn, AXLexpr, MerTKexpr, Tyro3expr)

    tps = range(0.0, 1440, length = 200)
    data = calcData(chn, AXLexpr, MerTKexpr, Tyro3expr, tps)

    pYAXL = TAMode.pY .* TAMode.recpSpecific[1]
    surfAXL = TAMode.surface .* TAMode.recpSpecific[1]
    totAXL = TAMode.total .* TAMode.recpSpecific[1]

    pYMerTK = TAMode.pY .* TAMode.recpSpecific[2]
    surfMerTK = TAMode.surface .* TAMode.recpSpecific[2]
    totMerTK = TAMode.total .* TAMode.recpSpecific[2]

    pYTyro3 = TAMode.pY .* TAMode.recpSpecific[3]
    surfTyro3 = TAMode.surface .* TAMode.recpSpecific[3]
    totTyro3 = TAMode.total .* TAMode.recpSpecific[3]

    pYA = (data * pYAXL)
    surfA = (data * surfAXL)
    totalA = (data * totAXL)

    pYM = (data * pYMerTK)
    surfM = (data * surfMerTK)
    totalM = (data * totMerTK)

    pYT = (data * pYTyro3)
    surfT = (data * surfTyro3)
    totalT = (data * totTyro3)

    A = hcat(pYA, surfA, totalA)
    M = hcat(pYM, surfM, totalM)
    T = hcat(pYT, surfT, totalT)
    plot(
        tps,
        [A, M, T],
        lw = 3,
        xlabel = "Time (min)",
        label = ["AXL" "AXL" "AXL" "MerTK" "MerTK" "MerTK" "Tyro3" "Tyro3" "Tyro3"],
        title = ["pY" "surf" "total"],
        layout = (1, 3),
        size = (950, 400),
    )
end

function calcData(chn, AXLexpr, MerTKexpr, Tyro3expr, tps)
    index = 201
    Ig2rev = get(chn, :Ig2rev)[1]

    x = get(chn, [:internalize, :pYinternalize, :sortF, :kRec, :kDeg, :xFwd, :gasCur, :AXLexpr])
    samp_params = hcat(x.internalize, x.pYinternalize, x.sortF, x.kRec, x.kDeg, x.xFwd, x.gasCur)
    params = vcat(samp_params[index, :], AXLexpr, MerTKexpr, Tyro3expr, Ig2rev[index], [1.0, 1.0, 1.8, 100.0])

    data = TAMode.runTAM(tps, params, 10)
    return data
end

"Plot Comp Model."
function plotComp(pp, tps)
    r = 1:100
    pY = TAMode.compTAM(tps, pp)
    cplot = Array{Float64}(undef, length(tps), length(tps), 100)

    for rr = 1:100
        for t = 1:length(tps)
            pYdata = view(pY, t, :, rr)
            cplot[:, rr] .= dot(pYdata, TAMode.pYc)
        end
    end

    plotpY = view(cplot, 1, 1, :)
    plot(r, plotpY, title = "Total Receptor pY", label = tps[1], lw = 3)

    if length(tps) > 1
        for tt = 2:length(tps)
            plotpY = view(cplot, tt, tt, :)
            plot!(r, plotpY, label = tps[tt], lw = 3)
        end
    end
    xlabel!("Radius")
    ylabel!("pY")

end

"Plot Comp Model Dimers."
function plotCompDimers(pp, tps, dimerVector)
    r = 1:100
    pY = TAMode.compTAM(tps, pp)
    dimerData = Array{Float64}(undef, length(tps), 100)

    for rr = 1:100
        for t = 1:length(tps)
            pYdata = view(pY, t, :, rr)
            dimerData[:, rr] .= dot(pYdata, dimerVector)
        end
    end

    plot_dimers = dimerData[1, :]
    plot(r, plot_dimers, label = tps[1], lw = 3)

    if length(tps) > 1
        for tt = 2:length(tps)
            plot_dimers = dimerData[tt, :]
            plot!(r, plot_dimers, label = tps[tt], lw = 3)
        end
    end

    title!("Dimer Formation")
    xlabel!("Radius")
    ylabel!("Dimer Concentration")

end

"Plot receptor-specific pY."
function plotRecpSpecific(pp, tps, recpVector)
    r = 1:100
    pY = TAMode.compTAM(tps, pp)
    pYrecp = TAMode.pYc .* recpVector

    cplot_recp = Array{Float64}(undef, length(tps), length(tps), 100)

    for rr = 1:100
        for t = 1:length(tps)
            pYdata = view(pY, t, :, r)
            cplot_recp[t, :, rr] .= dot(pYdata, pYrecp)
        end
    end

    plotpY_recp = view(cplot_recp, 1, 1, :)
    plot(r, plotpY_recp, title = "Receptor-Specific pY", label = tps[1], lw = 3)

    if length(tps) > 1
        for tt = 2:length(tps)
            plotpY_recp = view(cplot_recp, tt, tt, :)
            plot!(r, plotpY_recp, label = tps[tt], lw = 3)
        end
    end
    xlabel!("Radius")
    ylabel!("pY")

end


text:
using juliaAllotter
using Test

@testset "juliaAllotter.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: saurabhkm/juliaAllotter

max_issues_repo_name: saurabhkm/juliaAllotter

max_forks_repo_name: saurabhkm/juliaAllotter


text:
using RecyclePermanent
using MonteCarloMarkovKernels: simulateChainProgress, makeAMKernel, estimateBM,
  kde
using StaticArrays
using StatsBase
using LinearAlgebra
using JLD2
using Plots
Plots.gr()

include("../models/gaussianModel.jl")

@load "data/gaussianData.jld" ys θ
ys = ys[1:10]

d = 2

@inline function toGaussianθ(v::SVector{d, Float64})
  return Gaussianθ(v[1], v[2], θ.C, θ.D)
end

θ0 = SVector{d, Float64}(θ.A, θ.B)

const sigmaProp = SMatrix{d, d, Float64}(Matrix{Float64}(I, d, d))

@inline function lglogprior(theta::Gaussianθ)
  if theta.A < 0 || theta.A > 5 return -Inf end
  if theta.B < 0 || theta.B > 5 return -Inf end
  return 0.0
end

function makeSimpleLTD(ys::Vector{Float64}, N::Int64)
  rpio::RPIO{Float64Particle} = RPIO{Float64Particle}(N, length(ys), 1)
  function ltd(in::SVector{d, Float64})
    theta::Gaussianθ = toGaussianθ(in)
    lp::Float64 = lglogprior(theta)
    if lp == -Inf return -Inf end
    model::Model = makeGaussianModel(theta, ys)
    return simpleEstimate(model, rpio)
  end
end

function makeRecycleLTD(ys::Vector{Float64}, N::Int64)
  rpio::RPIO{Float64Particle} = RPIO{Float64Particle}(N, length(ys), 1)
  function ltd(in::SVector{d, Float64})
    theta::Gaussianθ = toGaussianθ(in)
    lp::Float64 = lglogprior(theta)
    if lp == -Inf return -Inf end
    model::Model = makeGaussianModel(theta, ys)
    return rpEstimate(model, rpio)
  end
end

function makeExactLTD(ys::Vector{Float64})
  function ltd(in::SVector{d, Float64})
    theta::Gaussianθ = toGaussianθ(in)
    lp::Float64 = lglogprior(theta)
    if lp == -Inf return -Inf end
    return gaussianExactLL(theta, ys)
  end
end

N = 40 * length(ys)
chainLength = 1000000

logtargetSimple = makeSimpleLTD(ys, N)
logtargetRecycle = makeRecycleLTD(ys, N)
logtargetExact = makeExactLTD(ys)

PSimple = makeAMKernel(logtargetSimple, sigmaProp)
PRecycle = makeAMKernel(logtargetRecycle, sigmaProp)
PExact = makeAMKernel(logtargetExact, sigmaProp)

Random.seed!(12345)

@time chainSimple = simulateChainProgress(PSimple, θ0, chainLength)
sar = PSimple(:acceptanceRate)

@time chainRecycle = simulateChainProgress(PRecycle, θ0, chainLength)
rar = PRecycle(:acceptanceRate)

@time chainExact = simulateChainProgress(PExact, θ0, chainLength)
ear = PExact(:acceptanceRate)

vsSimple = (i->(x->x[i]).(chainSimple)).(1:d)
vsRecycle = (i->(x->x[i]).(chainRecycle)).(1:d)
vsExact = (i->(x->x[i]).(chainExact)).(1:d)

plot(kde(vsSimple[1], sar))
plot!(kde(vsExact[1], ear))
plot!(kde(vsRecycle[1], rar))

plot(kde(vsSimple[2], sar))
plot!(kde(vsExact[2], ear))
plot!(kde(vsRecycle[2], rar))

contour(kde(vsSimple[1], vsSimple[2], sar))
contour!(kde(vsExact[1], vsExact[2], ear))
contour!(kde(vsRecycle[1], vsRecycle[2], rar))

plot(autocor(vsSimple[1]))
plot!(autocor(vsExact[1]))
plot!(autocor(vsRecycle[1]))

plot(autocor(vsSimple[2]))
plot!(autocor(vsExact[2]))
plot!(autocor(vsRecycle[2]))

avarsSimple = estimateBM.(vsSimple)
avarsRecycle = estimateBM.(vsRecycle)
avarsExact = estimateBM.(vsExact)