text:
mutable struct abstract_var
end

mutable struct context
end

mutable struct empinfo
end

mutable struct equtree
end

mutable struct equnode
end

mutable struct jamsd_var
end

mutable struct jamsd_equ
end

mutable struct jamsd_sp_matrix
end

mutable struct mathprgm
end

mutable struct ovf_definition
end

mutable struct jamsd_options
end

mutable struct equil
end

mutable struct option
	name::Vector{Cchar}
	typ::Cint
	value::Cdouble
end

# TODO with latest Julia
if VERSION < v"0.7"
	iswin = is_windows()
	Cvoid = Void
else
	iswin = Sys.iswindows()
end

#const libjamsd = iswin ? "jamsd" : "libjamsd"
const term_str = iswin ? "\r\n" : "\n"

function ctx_add_lin_var(ctx::Ptr{context}, eidx, vidx, coeff::Cdouble)
	equ = ctx_getequ(ctx, eidx)
	res = ccall((:equ_add_var, libjamsd), Cint, (Ptr{context}, Ptr{jamsd_equ}, Cint, Cdouble), ctx, equ, vidx, coeff)
	res != 0 && error("return code $res from JAMSD")
end

function ctx_create(n, m)
	ctx = ccall((:ctx_alloc, libjamsd), Ptr{context}, (Cuint,), 2)
	res = ccall((:model_reserve_eqns, libjamsd), Cint, (Ptr{context}, Ptr{Cvoid}, Cuint), ctx, C_NULL, m);
	res != 0 && error("return code $res from JAMSD")
	res = ccall((:model_reserve_vars, libjamsd), Cint, (Ptr{context}, Ptr{Cvoid}, Cuint), ctx, C_NULL, n);
	res != 0 && error("return code $res from JAMSD")
	res = ccall((:ctx_resize, libjamsd), Cint, (Ptr{context}, Cuint, Cuint), ctx, n, m)
	res != 0 && error("return code $res from JAMSD")
	return ctx
end

function ctx_getvar(ctx, idx)
	return ccall((:ctx_getvar, libjamsd), Ptr{jamsd_var}, (Ptr{context}, Cint), ctx, idx)
end

function ctx_m(ctx)
	return ccall((:ctx_m, libjamsd), Cint, (Ptr{context},), ctx)
end

function ctx_n(ctx)
	return ccall((:ctx_n, libjamsd), Cint, (Ptr{context},), ctx)
end

function ctx_setvarnames(ctx, names::Vector{String})
	res = ccall((:myo_set_varnames, libjamsd), Cint, (Ptr{context}, Ptr{Ptr{Cchar}}, Cuint), ctx, names, length(names))
	res != 0 && error("return code $res from JAMSD")
end

function hack_last_vidx(ctx)
	return ccall((:model_total_n, libjamsd), Csize_t, (Ptr{context},), ctx) - 1
end

function hack_exportempinfo(ctx, ctx_mtr, emp)
	res = ccall((:hack_exportempinfo, libjamsd), Cint, (Ptr{context}, Ptr{context}, Ptr{empinfo}), ctx, ctx_mtr, emp)
	res != 0 && error("return code $res from JAMSD")
end

function hack_solver_log()
	ccall((:hack_solver_log, libjamsd), Cvoid, ())
end

function ctx_setvarlone(ctx::Ptr{context}, idx, val::Cdouble)
	res = ccall((:ctx_setvarlone, libjamsd), Cint, (Ptr{context}, Cint, Cdouble), ctx, idx, val)
	res != 0 && error("return code $res from JAMSD")
end

function ctx_getvarval(ctx::Ptr{context}, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getvarlone, libjamsd), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from JAMSD")
	return val.x
end

function ctx_getvarmult(ctx::Ptr{context}, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getvarmone, libjamsd), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from JAMSD")
	return val.x
end

function ctx_getequ(ctx, idx)
	equ = ccall((:ctx_getequ, libjamsd), Ptr{jamsd_equ}, (Ptr{context}, Cint), ctx, idx)
end

function ctx_getequval(ctx, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getequlone, libjamsd), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from JAMSD")
	return val.x
end

function ctx_getmultiplierval(ctx::Ptr{context}, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getequmone, libjamsd), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from JAMSD")
	return val.x
end

function emp_create(ctx)
	return ccall((:empinfo_alloc, libjamsd), Ptr{empinfo}, (Ptr{context},), ctx)
end

function emp_create_equil(max_mp)
	return ccall((:mp_equil_alloc, libjamsd), Ptr{equil}, (Cuint,), max_mp)
end

function emp_delete(emp::Ptr{empinfo})
	return ccall((:hack_empinfo_dealloc, libjamsd), Cvoid, (Ptr{empinfo},), emp)
end

function emp_add_mp_mp(mp_parent, mp)
	res = ccall((:mathprgm_addmp, libjamsd), Cint, (Ptr{mathprgm}, Ptr{mathprgm}), mp_parent, mp)
	res != 0 && error("return code $res from JAMSD")
end

function emp_add_mp_equil(mp_parent, mpe)
	res = ccall((:mathprgm_addequil, libjamsd), Cint, (Ptr{mathprgm}, Ptr{equil}), mp_parent, mpe)
	res != 0 && error("return code $res from JAMSD")
end

function emp_add_equil(emp)
	mpe = Ref{Ptr{equil}}(C_NULL)
	res = ccall((:empinfo_add_equil, libjamsd), Cint, (Ptr{empinfo}, Ref{Ptr{equil}}), emp, mpe)
	res != 0 && error("return code $res from JAMSD")

	return mpe.x
end

function emp_equil_add(mpe, mp)
	res = ccall((:mp_equil_add, libjamsd), Cint, (Ptr{equil}, Ptr{mathprgm}), mpe, mp)
	res != 0 && error("return code $res from JAMSD")
end

function emp_hack(emp)
	res = ccall((:hack_ag_addfinish, libjamsd), Cint, (Ptr{empinfo},), emp)
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_ensure(emp, nb)
	res = ccall((:empinfo_ensure, libjamsd), Cint, (Ptr{empinfo}, Cuint), emp, nb)
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_alloc(emp, ctx)
	mp = ccall((:mathprgm_alloc, libjamsd), Ptr{mathprgm}, (Ptr{empinfo}, Ptr{context}), emp, ctx)
	return mp
end

function emp_mp_start(mp, typ)
	res = ccall((:mathprgm_addstart, libjamsd), Cint, (Ptr{mathprgm}, Cuint), mp, typ)
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_objdir(mp, sense)
	res = ccall((:mathprgm_addobjdir, libjamsd), Cint, (Ptr{mathprgm}, Cint), mp, sense_to_jamsd[sense])
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_objequ(mp, idx)
	res = ccall((:mathprgm_addobjequ, libjamsd), Cint, (Ptr{mathprgm}, Cint), mp, idx)
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_objvar(mp, idx)
	res = ccall((:mathprgm_addobjvar, libjamsd), Cint, (Ptr{mathprgm}, Cint), mp, idx)
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_var(mp, idx)
	res = ccall((:mathprgm_addvar, libjamsd), Cint, (Ptr{mathprgm}, Cint), mp, idx)
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_constraint(mp, idx)
	res = ccall((:mathprgm_addconstraint, libjamsd), Cint, (Ptr{mathprgm}, Cint), mp, idx)
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_vipair(mp, eidx, vidx)
	res = ccall((:mathprgm_addvipair, libjamsd), Cint, (Ptr{mathprgm}, Cint, Cint), mp, eidx, vidx)
	res != 0 && error("return code $res from JAMSD")
end

function emp_mp_getobjvar(mp)
	return ccall((:mathprgm_getobjvar, libjamsd), Cint, (Ptr{mathprgm},), mp)
end

function emp_mp_getobjdir(mp)
	return ccall((:mathprgm_getobjvar, libjamsd), Cint, (Ptr{mathprgm},), mp)
end

function emp_mp_print(mp, ctx)
	ccall((:mathprgm_print, libjamsd), Cint, (Ptr{mathprgm}, Ptr{context}), mp, ctx)
end
#function emp_mp_to_agent(ctx, emp)
#	res = ccall((:mathprgm_to_agent, libjamsd), Cint, (Ptr{context}, Ptr{empinfo}), ctx, emp)
#	res != 0 && error("return code $res from JAMSD")
#end

function emp_report_values(emp, ctx)
	res = ccall((:empinfo_report_values, libjamsd), Cint, (Ptr{empinfo}, Ptr{context}), emp, ctx)
	res != 0 && error("return code $res from JAMSD")
end

function emp_set_root(emp, mpe::Ptr{equil})
	res = ccall((:empinfo_set_emproot_mpe, libjamsd), Cint, (Ptr{empinfo}, Ptr{equil}), emp, mpe)
	res != 0 && error("return code $res from JAMSD")
end

function emp_set_root(emp, mp::Ptr{mathprgm})
	res = ccall((:empinfo_set_emproot_mp, libjamsd), Cint, (Ptr{empinfo}, Ptr{mathprgm}), emp, mp)
	res != 0 && error("return code $res from JAMSD")
end

function emp_ovf(emp, name, ovf_vidx, args)
	ovf_def = Ref{Ptr{ovf_definition}}(C_NULL)
	res = ccall((:ovf_add, libjamsd), Cint, (Ptr{empinfo}, Cstring, Cint, Cuint, Ptr{Cint}, Ref{Ptr{ovf_definition}}),
							                               emp, name, ovf_vidx, length(args), args, ovf_def)
	res != 0 && error("return code $res from JAMSD")

	return ovf_def.x
end

function emp_ovf_param(ovf_def, param_name, scalar::Number)
	res = ccall((:ovf_param_add_scalar, libjamsd), Cint, (Ptr{ovf_definition}, Cstring, Cdouble), ovf_def, param_name, scalar)
	res != 0 && error("return code $res from JAMSD")
end

function emp_ovf_param(ovf_def, param_name, arr::Vector)
  arrC = Vector{Cdouble}(arr)
	res = ccall((:ovf_param_add_vector, libjamsd), Cint, (Ptr{ovf_definition}, Cstring, Cuint, Ptr{Cdouble}), ovf_def, param_name, length(arr), arrC)
	res != 0 && error("return code $res from JAMSD")
end

function emp_ovf_check(ovf_def)
	res = ccall((:ovf_check, libjamsd), Cint, (Ptr{ovf_definition},), ovf_def)
end

function equtree_var(ctx, tree, node, idx, coeff)
	res = ccall((:equtree_var, libjamsd), Cint, (Ptr{context}, Ptr{equtree}, Ref{Ref{Ptr{equnode}}}, Cint, Cdouble),
		ctx,
		tree,
		node,
		idx,
		coeff)
	res != 0 && error("return code $res from JAMSD")
end

function equtree_cst(ctx, tree, node, value)
	res = ccall((:equtree_cst, libjamsd), Cint, (Ptr{context}, Ptr{equtree}, Ref{Ref{Ptr{equnode}}}, Cdouble),
		ctx,
		tree,
		node,
		value)
	res != 0 && error("return code $res from JAMSD")
end

function equtree_arithm(tree, node, opcode, nb)
	res = ccall((:equtree_arithm, libjamsd), Cint, (Ptr{equtree}, Ref{Ref{Ptr{equnode}}}, Cuint, Cuint),
		tree,
		node,
		opcode,
		nb)
	res != 0 && error("return code $res from JAMSD")
end

function equtree_call(ctx, tree, node, fndata)
	res = ccall((:equtree_call, libjamsd), Cint, (Ptr{context}, Ptr{equtree}, Ref{Ref{Ptr{equnode}}}, Cuint, Cuint),
		ctx,
		tree,
		node,
		fndata[1],
		fndata[2]
		)
	res != 0 && error("return code $res from JAMSD")
end

function equtree_get_root_addr(tree::Ptr{equtree}, node::Ref{Ref{Ptr{equnode}}})
	res = ccall((:equtree_get_root_addr, libjamsd), Cint, (Ptr{equtree}, Ref{Ref{Ptr{equnode}}}),
		tree,
		node)
	res != 0 && error("return code $res from JAMSD")
end

function equtree_umin(ctx, tree, node)
	res = ccall((:equtree_umin, libjamsd), Cint, (Ptr{equtree}, Ref{Ref{Ptr{equnode}}}),
		tree,
		node)
	res != 0 && error("return code $res from JAMSD")
end

function equnode_get_child_addr(node::Ptr{equnode}, i::Int)
	child = Ref{Ref{Ptr{equnode}}}(C_NULL)
	res = ccall((:equnode_get_child_addr, libjamsd), Cint, (Ptr{equnode}, Ref{Ref{Ptr{equnode}}}, Cuint),
		node,
		child,
		i)
	res != 0 && error("return code $res from JAMSD")
	return child
end

function equnode_deref(node::Ref{Ref{Ptr{equnode}}})
	return ccall((:p2deref, libjamsd), Ptr{equnode}, (Ref{Ref{Ptr{equnode}}},), node)
end

function ctx_get_solvername(ctx::Ptr{context})
	str = "?"^256
	res = ccall((:ctx_getsolverstr, libjamsd), Cint, (Ptr{context}, Cstring, Cuint), ctx, str, length(str))
	res != 0 && error("return code $res from JAMSD")
	return str
end

function jamsd_add_box_var(ctx::Ptr{context}, lower::Cdouble, upper::Cdouble)
	res = ccall((:model_add_box_var, libjamsd), Cint, (Ptr{context}, Cdouble, Cdouble), ctx, lower, upper)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_add_free_var(ctx::Ptr{context}, nb)
	res = ccall((:model_add_free_vars, libjamsd), Cint, (Ptr{context}, Cuint), ctx, nb)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_add_neg_var(ctx::Ptr{context}, nb)
	res = ccall((:model_add_neg_vars, libjamsd), Cint, (Ptr{context}, Cuint), ctx, nb)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_add_pos_var(ctx::Ptr{context}, nb)
	res = ccall((:model_add_pos_vars, libjamsd), Cint, (Ptr{context}, Cuint), ctx, nb)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_avar(size, indices)
	indicesC = Vector{Cint}(indices)
	return ccall((:avar_alloc, libjamsd), Ptr{abstract_var}, (Cuint, Ptr{Cint}), size, indicesC)
end

function jamsd_avar_free(avar::Ptr{abstract_var})
	return ccall((:avar_free, libjamsd), Cvoid, (Ptr{abstract_var},), avar)
end

function jamsd_decl_eqn(ctx::Ptr{context}, idx)
	minn = Ref{Cint}(-1)
	res = ccall((:model_add_eqn_empty, libjamsd), Cint, (Ptr{context}, Ref{Cint}, Ptr{Cvoid}, Cint), ctx, minn, C_NULL, 0)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_equ_add_lin_tree(ctx::Ptr{context}, eidx, qvals, avar::Ptr{abstract_var}, coeff)
	equ = ctx_getequ(ctx, eidx)
	res = ccall((:equ_add_lin_tree, libjamsd), Cint, (Ptr{context}, Ptr{jamsd_equ}, Ref{Cdouble}, Ptr{abstract_var}, Cdouble),
							ctx, equ, qvals, avar, coeff)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_equ_add_quadratic(ctx::Ptr{context}, eidx, mat::Ptr{jamsd_sp_matrix}, avar::Ptr{abstract_var}, coeff)
	equ = ctx_getequ(ctx, eidx)
	res = ccall((:equ_add_quadratic, libjamsd), Cint, (Ptr{context}, Ptr{jamsd_equ}, Ptr{jamsd_sp_matrix}, Ptr{abstract_var}, Cdouble),
							ctx, equ, mat, avar, coeff)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_mat_coo(ridx, cidx, vals)
	# beaware of dragons! --xhub
	ridxC = Vector{Cint}(ridx)
	cidxC = Vector{Cint}(cidx)
	# m and n are not really needed? --xhub
	mat = ccall((:empmat_triplet, libjamsd), Ptr{jamsd_sp_matrix}, (Cuint, Cuint, Cuint, Ptr{Cint}, Ptr{Cint}, Ref{Cdouble}),
							0, 0, length(vals), ridxC, cidxC, vals)
	return mat
end

function jamsd_mat_free(mat)
	ccall((:empmat_free, libjamsd), Cvoid, (Ptr{jamsd_sp_matrix},), mat)
end

function jamsd_set_objeqn(ctx::Ptr{context}, idx)
	res = ccall((:model_setobjequ, libjamsd), Cint, (Ptr{context}, Cint), ctx, idx)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_set_rhs(ctx::Ptr{context}, idx, val)
	res = ccall((:ctx_setrhs, libjamsd), Cint, (Ptr{context}, Cint, Cdouble), ctx, idx, val)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_set_equtype(ctx::Ptr{context}, idx, rel)
	res = ccall((:ctx_setequtype, libjamsd), Cint, (Ptr{context}, Cint, Cuint), ctx, idx, rel)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_set_vartype(ctx::Ptr{context}, idx, typ)
	res = ccall((:ctx_setvartype, libjamsd), Cint, (Ptr{context}, Cint, Cuint), ctx, idx, typ)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_get_treedata(ctx, i::Int)
	tree = ccall((:myo_getequtree, libjamsd), Ptr{equtree}, (Ptr{context}, Cint), ctx, i)
	node = Ref{Ref{Ptr{equnode}}}(C_NULL)
	res = equtree_get_root_addr(tree, node)
	res != 0 && error("return code $res from JAMSD")
	return (tree, node)
end

function jamsd_options_alloc()
	return ccall((:hack_options_alloc, libjamsd), Ptr{jamsd_options}, ())
end

function jamsd_options_dealloc(o::Ptr{jamsd_options})
	ccall((:hack_options_dealloc, libjamsd), Cvoid, (Ptr{jamsd_options},), o)
end

function jamsd_option_set(opt::Ptr{jamsd_options}, k::String, v::Bool)
	res = ccall((:option_set_b, libjamsd), Cint, (Ptr{jamsd_options}, Cstring, Cint), opt, k, v)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_option_set(opt::Ptr{jamsd_options}, k::String, v::Int)
	res = ccall((:option_set_i, libjamsd), Cint, (Ptr{jamsd_options}, Cstring, Cint), opt, k, v)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_option_set(opt::Ptr{jamsd_options}, k::String, v::Cdouble)
	res = ccall((:option_set_d, libjamsd), Cint, (Ptr{jamsd_options}, Cstring, Cdouble), opt, k, v)
	res != 0 && error("return code $res from JAMSD")
end

function jamsd_option_set(opt::Ptr{jamsd_options}, k::String, v::String)
	res = ccall((:option_set_s, libjamsd), Cint, (Ptr{jamsd_options}, Cstring, Cstring), opt, k, v)
	res != 0 && error("return code $res from JAMSD")
end

function ctx_dealloc(o::Ptr{context})
	ccall((:hack_ctx_dealloc,  libjamsd), Cvoid, (Ptr{context},), o)
end


text:
# Allow one and two-argument forms
cost(::EnergyModel, d::Device) = cost(d)
cost(::Device) = error("Not implemented")

expression(::Expression, m::EnergyModel, d::Device) = error("Not implemented")
expression(::Cost, m::EnergyModel, d::Device) = cost(d)
function expression(::NodalActivePower, m::EnergyModel, b::Bus)
    T = axis(m, :snapshots)
    buses = axis(m, b)
    p = AxisArray(zeros(AffExpr, length(buses), length(T)), buses, T)
    for d in devices(model(b))
        addtoexpr!(p, NodalActivePower(), m, d)
    end
    p
end

"""
    addto!(jm::JuMP.AbstractModel, m::EnergyModel)
    addto!(jm::JuMP.AbstractModel, m::EnergyModel, d::Device)

Add the variables and constraints defining the energy model `m` to the jump
model `jm` and set the objective function. If an device `d` is specified, the
device is added alone.
"""
function addto! end

# Wraps the AbstractModel in a ModelView to store variables and constraints in m.objects[<class>]
function addto!(jm::JuMP.AbstractModel, m::EnergyModel, c::Component)
    addto!(ModelView(jm, c), m, c)
end
function addto!(jm::JuMP.AbstractModel, m::EnergyModel, d::Device)
    addto!(ModelView(jm, d), m, demote_formulation(m, d))
end

include("energymodel.jl")
include("subnetwork.jl")
include("bus.jl")

## OnePort

"Connected to exactly one `Bus`, determined by :bus attribute"
abstract type OnePort{DF<:DeviceFormulation} <: Device{DF} end
busattributes(d::OnePort) = (:bus,)

## Defaults for OnePort
cost(d::OnePort) = sum(d[:marginal_cost] .* d[:p]) + sum(d[:capital_cost] .* (d[:p_nom] .- getparam(d, :p_nom)))
function addtoexpr!(injection::AxisArray, ::NodalActivePower,
                    m::EnergyModel, d::OnePort)
    B, T = AxisArrays.axes(injection)
    D = axis(m, d)

    buses = get(d, :bus, D)
    p = get(d, :p, D, T)

    for (i, idx) in enumerate(indexin(buses, B)), t in T
        !isnothing(idx) && add_to_expression!(injection[idx,t], p[i,t])
    end
end

include("generator.jl")
include("storageunit.jl")
include("store.jl")
include("load.jl")

# Branch
"Connected to exactly two `Bus`es, determined by :bus0 and :bus1 attribute"
abstract type Branch{DF<:DeviceFormulation} <: Device{DF} end

# ActiveBranch
abstract type ActiveBranch{DF<:DeviceFormulation} <: Branch{DF} end

busattributes(d::Branch) = (:bus0, :bus1)
function p(d::Branch)
    p = d[:p]
    ((b,t)->p[b,t],   # :bus0
     (b,t)->-p[b,t])  # :bus1
end

function addtoexpr!(injection::AxisArray, ::NodalActivePower,
                    m::EnergyModel, d::Branch)
    B, T = AxisArrays.axes(injection)
    L = axis(m, d)

    bus0 = get(d, :bus0, L)
    bus1 = get(d, :bus1, L)
    p = get(d, :p, L, T)

    for (i, idx0, idx1) in zip(1:length(L), indexin(bus0, B), indexin(bus1, B)), t in T
        !isnothing(idx0) && add_to_expression!(injection[idx0,t], p[i,t])
        !isnothing(idx1) && add_to_expression!(injection[idx1,t], -1, p[i,t])
    end
end

include("link.jl")
include("passivebranch.jl")

ext: jl

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
module Metrics
# Module for DL Metrics

export mae, mse, msle, male, r2_score, adjusted_r2_score

export binary_accuracy, confusion_matrix, categorical_accuracy, sparse_categorical, top_k_categorical, top_k_sparse_categorical, precision, recall, sensitivity, detection_rate, f_beta_score, specificity, false_alarm_rate, cohen_kappa, statsfromTFPN, classwise_stats, global_stats

export IoU, PSNR

export bleu_score, rouge, rouge_l_summary_level

export ranking_stats_k, avg_precision

export report_stats

using StatsBase, DataFrames
using DataStructures: OrderedDict

include("./CV_Metrics/CVMetrics.jl")
include("./NLP_Metrics/NLPMetrics.jl")
include("Classification.jl")
include("Regression.jl")
include("Ranking_n_Statistical.jl")
include("utils.jl")

end


text:
n_samples = 100_000

for d in [MultivariateNormal([0.0, 0.0], [1.0 0.9; 0.9 1.0]),
          Dirichlet([2.0, 1.0]),
          Multinomial(1, [0.5, 0.5])]
    # Test rand
    s = rand(d)
    s = rand(d, 10)
    A = Array(eltype(s), 2, n_samples)
    rand!(d, A)

    # Test mean under CLT
    mu_hat = mean(A, 2)
    @assert norm(mu_hat - mean(d), Inf) < 0.1

    # Test cov under CLT
    sigma_hat = cov(A')
    @assert norm(sigma_hat - var(d), Inf) < 0.1

    # Test entropy and logpdf using AEP
    en_hat = -mean(logpdf(d, A))
    en = entropy(d)
    @assert abs(en_hat - en) < 0.1

    # Test entropy and pdf using AEP
    en_hat = -mean(log(pdf(d, A)))
    en = entropy(d)
    @assert abs(en_hat - en) < 0.1
end


text:
using FGES
using CSV, DataFrames
using Combinatorics
using Statistics

#Generate a dataset
numFeatures = 100
numObservations = 10000
data = zeros(numObservations, numFeatures)

#number of previous features to use in caluclating the next feature
n=5

for i in 1:numFeatures
    if i ≤ n
        data[:,i] = randn(numObservations)
    else
        data[:,i] = sum(rand()*data[:,i-j] for j∈1:n) + randn(numObservations)
        #Stop the data from exploding
        data[:,i] ./= mean(data[:,i])
    end
end

#Run the julia version of FGES and time it
gJulia = @time fges(data)

#Save the dataset to run for java
df = DataFrame(data, :auto)
CSV.write("test/testDatasets/features$(numFeatures)_observations$(numObservations).txt", df ,delim=" ")

#Read and parse the java results
gJava = PDAG(numFeatures)
javaFilepath = "test/testDatasets/features$(numFeatures)_observations$(numObservations)_java.txt"

javaResult = filter(line -> occursin(r"x\d+ \D{3} x\d+",line),readlines(open(javaFilepath)))
javaResult = split.(javaResult," ")

javaEdges = [parse.(Int, filter.(isdigit, line[[2,4]])) for line in javaResult]

for (i,edge) in enumerate(javaEdges)
    addedge!(gJava,edge...,directed=javaResult[i][3]=="-->" ? true : false) #not worried about orientation for now
end

#Generate the true graph
gTrue = PDAG(numFeatures)

for i=n+1:numFeatures
    for j=1:n
        addedge!(gTrue,i-j, i)
    end
end

#Calculate precision and recall
function contingencyTable(g,gTrue)

    continTable = zeros(Int,2,2)

    for (x,y) in combinations(vertices(gTrue),2)
        trueEdge = isadjacent(gTrue,x,y)
        estEdge = isadjacent(g,x,y)

        if trueEdge & estEdge #true positive
            continTable[1,1] += 1
        elseif !trueEdge & !estEdge #true negative
            continTable[2,2] += 1
        elseif !trueEdge & estEdge #false positive
            continTable[2,1] += 1
        else #false negative
            continTable[1,2] += 1
        end
    end

    return continTable
end

modelPrecision(continTable) = continTable[1,1] / sum(continTable[:,1])
modelRecall(continTable) = continTable[1,1] / sum(continTable[1,:])

continTableJulia = contingencyTable(gJulia,gTrue)
modelPrecision(continTableJulia)
modelRecall(continTableJulia)

continTableJava = contingencyTable(gJava,gTrue)
modelPrecision(continTableJava)
modelRecall(continTableJava)

#Calculate precision and recall
function contingencyTableOrient(g,gTrue)

    continTable = zeros(Int,2,2)

    for (x,y) in permutations(vertices(gTrue),2)
        # trueAdj = isadjacent(gTrue,x,y)
        # estAdj = isadjacent(g,x,y)

        # trueDir = isparent(gTrue,x,y)
        # estDir = isparent(g,x,y)
        trueDir = Edge(x,y,true) ∈ edges(gTrue)
        estDir = Edge(x,y,true) ∈ edges(g)

        # if (trueAdj & estAdj) & (!trueUndir & !estUndir) #true positive
        #     continTable[1,1] += 1
        # elseif (!trueAdj & !estAdj) | (trueUndir & estUndir) #true negative
        #     continTable[2,2] += 1
        # elseif (!trueAdj | trueUndir) & (estAdj & !estUndir) #false positive
        #     continTable[2,1] += 1
        # else #false negative
        #     continTable[1,2] += 1
        # end

        if trueDir & estDir
            continTable[1,1] += 1
        elseif !trueDir & !estDir
            continTable[2,2] += 1
        elseif !trueDir & estDir
            continTable[2,1] += 1
        else
            continTable[1,2] += 1
        end
    end

    return continTable
end

continTableJulia = contingencyTableOrient(gJulia,gTrue)
modelPrecision(continTableJulia)
modelRecall(continTableJulia)

continTableJava = contingencyTableOrient(gJava,gTrue)
modelPrecision(continTableJava)
modelRecall(continTableJava)


text:
using Rotolo

@session test style=>"color:red"
@redirect String
"red"

@container red class=>"redbox" style=>"color:red"
"red"
@container red id=>"test" class=>"bluebox" style=>"color:blue"
"red"

@container blue style=>"color:green"
"green"
@container blue.red style=>"color:yellow"
"yellow"
@container blue.green style=>"color:blue"
"blue"

@container red3 class=>"red2" style=>"color:red"
"red"

theme = "font-size:large;" ;

@container red.large style => theme
"bigger"

@style "bigger" style => "font-size:x-large;font-style:italic"

@style 456.45 style=>"font-size:x-small"

@style "test" style=>"font-style:italic" class=>"test" id=>"youpi"


text:
const CPetscViewer = Ptr{Cvoid}

"""
    AbstractViewer{PetscLib <: PetscLibType}

Abstract type of PETSc viewer.

# External Links
$(_doc_external("Viewer/PetscViewer"))
"""
abstract type AbstractViewer{PetscLib <: PetscLibType} end

"""
    ViewerStdout(petsclib, comm = MPI.COMM_SELF)

Create an ASCII `PetscViewer` for the `comm`

# External Links
$(_doc_external("Viewer/PETSC_VIEWER_STDOUT_"))
"""
mutable struct ViewerStdout{PetscLib} <: AbstractViewer{PetscLib}
    ptr::CPetscViewer
end

@for_petsc function ViewerStdout(
    ::$UnionPetscLib,
    comm::MPI.Comm,
)
    ptr = ccall(
        (:PETSC_VIEWER_STDOUT_, $petsc_library),
        CPetscViewer,
        (MPI.MPI_Comm,),
        comm,
    )
    return ViewerStdout{$PetscLib}(ptr)
end

@for_petsc function Base.push!(
    viewer::AbstractViewer{$PetscLib},
    format::PetscViewerFormat,
)
    @chk ccall(
        (:PetscViewerPushFormat, $petsc_library),
        PetscErrorCode,
        (CPetscViewer, PetscViewerFormat),
        viewer,
        format,
    )
    return nothing
end

@for_petsc function Base.pop!(viewer::AbstractViewer{$PetscLib})
    @chk ccall(
        (:PetscViewerPopFormat, $petsc_library),
        PetscErrorCode,
        (CPetscViewer,),
        viewer,
    )
    return nothing
end

function with(f, viewer::AbstractViewer, format::PetscViewerFormat)
    push!(viewer, format)
    try
        f()
    finally
        pop!(viewer)
    end
end

# ideally we would capture the output directly, but this looks difficult
# easiest option is to redirect stdout
# based on suggestion from https://github.com/JuliaLang/julia/issues/32567
function _show(io::IO, obj)
    old_stdout = stdout
    try
        rd, = redirect_stdout()
        view(obj)

        # Since not all MPI ranks are guaranteed to print we put in a newline
        # that we remove with the write since readavailable will hang if there
        # is no data in the stream
        println()

        Libc.flush_cstdio()
        flush(stdout)
        write(io, readavailable(rd)[1:end-1])
    finally
        redirect_stdout(old_stdout)
    end
    return nothing
end

#=
# PETSc_jll isn't built with X support
mutable struct ViewerDraw <: AbstractViewer
    ptr::CPetscViewer
    comm::MPI.Comm
end
function ViewerDraw(comm::MPI.Comm)
    ptr = ccall((:PETSC_VIEWER_DRAW_, libpetsc), CPetscViewer, (MPI.MPI_Comm,), comm)
    return ViewerDraw(ptr, comm)
end
=#


text:
# Paralel network simulation

# Add process
using Distributed
nprocs() == Sys.CPU_THREADS || addprocs(Sys.CPU_THREADS - nprocs() + 1)

# Load packages
@everywhere begin
    using Pkg
    Pkg.activate(joinpath(Pkg.envdir(), "dev-env"))
end

@everywhere begin
    import Cairo, Fontconfig
    using GraphPlot
    using Graphs
    using SimpleWeightedGraphs
    using LinearAlgebra
    using DynamicalNetworks
    using DifferentialEquations
end

@everywhere begin
    σr = collect(0.0:0.1:1.0)
    function prob_func(prob, i, repeat)
        @show σr[i]
        prob.p[1].noisemat .*= σr[i]
        prob
    end

    # Settings
    nbits = 10
    tbit = 50.0
    bits = rand(Bool, nbits)
    ti = 0.0
    tf = nbits * tbit
    tswitch = ti:tbit:tf-tbit
    srcnode1, dstnode1 = 5, 8     # Control link nodes
    high = 100.0            # Control coupling strength
    low = eps()             # Control coupling strength
    condition(u, t, integrator) = integrator.t ∈ tswitch
    function action!(integrator)
        net, bits = integrator.sol.prob.p
        if popfirst!(bits)
            add_edge!(net.graph, srcnode1, dstnode1, high)
        else
            add_edge!(net.graph, srcnode1, dstnode1, low)
        end
    end
end

# Topology
n = 8
graph = cycle_graph(n) |> SimpleWeightedGraph
add_edge!(graph, 2, 7)
add_edge!(graph, 3, 6)
ϵ = 20.0

# Network
σ = σr[1]
H = [
    1 0 0 0
    -1 0 0 0
    0 1 0 0
    0 -1 0 0
    0 0 1 0
    0 0 -1 0
    0 0 0 1
    0 0 0 -1
] * σ
P = diagm([1.0, 0.0, 0.0])
net = SDENetwork(Lorenz, ϵ, graph, H, P)
callback = DiscreteCallback(condition, action!)

# Simulation
prob = netproblem(net, (ti, tf), copy(bits), noise_rate_prototype = get_noise_rate_prototype(net), callback = callback, tstops = tswitch)
ensembleprob = EnsembleProblem(prob, prob_func = prob_func)
sim = solve(ensembleprob, SRA(), EnsembleDistributed(), trajectories = length(σr))

# Plots
using PyPlot
pygui(true)
close("all")
getwave(x, i, d) = getindex.(x, (i - 1) * d + 1)
figure()
d = 3
subplot(221)
for (i, sol) in enumerate(sim[[10, 1]])
    plot(sol.t, abs.(getwave(sol.u, 1, d) - getwave(sol.u, 8, d)), label = "$i")
    axvline.(collect(ti:tbit:tf), color = "red", linestyle = "dashed")
    legend()
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 28042615-8f2b-41fa-9b14-28d69180debb
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 3ec6aae0-c233-11eb-355d-3dfe780e33bc
begin
	using DataFrames
	using JuMP
	using Ipopt
	using StatsPlots
	plotly()
end

# ╔═╡ a6997248-aee3-4ebb-8a1e-f791ad503aa8
begin
	n = 100
	μ1 = 10
	μ2 = 15

	df = append!(
		DataFrame(x1=randn(n).+μ1, x2=randn(n).+μ1, y=ones(n)),
		DataFrame(x1=randn(n).+μ2, x2=randn(n).+μ2, y=-ones(n))
	)

	p = @df df scatter(:x1, :x2, group=:y, legend=:topleft)
end

# ╔═╡ bf7d0459-4b47-4839-800d-01bdcd500633
md"
$arg\ min_{𝐰, b}\ \frac{1}{2}𝐰^T𝐰$

$\text{subject to }\forall i, y_i(𝐰^T𝐱_i+b) \geq 1$
"

# ╔═╡ 27dfce4e-3ccd-4221-8eca-144ae87e2998
md"
respect to

$arg\ min_{𝐰, b}\ \frac{1}{2}𝐮^T 𝐐 𝐮 + 𝐏^T 𝐮$

$\text{subject to }\forall i, 𝐚_i^T 𝐮) \geq c_i$

Therefore,

$𝐮 = \begin{bmatrix}b \\ 𝐰\end{bmatrix}$
$𝐐 = \begin{bmatrix}0 & 0^T_d \\ 0^T_d & 𝐈_d\end{bmatrix}$
$𝐏 = 0_{d+1}$
$𝐚^T_i = y_i\begin{bmatrix}1 & 𝐱^T_i\end{bmatrix}$
$c_i = 1$
"

# ╔═╡ 5056cc0a-f832-47b3-a6a5-f7ee2d7bc5ba
begin
	x = Matrix(df[:, 1:2])
	y = Vector(df[:, 3])

	function solve_svm(x, y)
		svm = Model(Ipopt.Optimizer)
		@variable(svm, w[1:2])
		@variable(svm, b)
		@objective(svm, Min, 0.5 * w' * w)
		@constraint(svm, y.*(x * w .+ b) .>= 1)
		optimize!(svm)

		@show objective_value(svm)

		return value.(w), value(b)
	end

	w, b = solve_svm(x, y)
end

# ╔═╡ e5334dab-56b2-4dd0-a12c-e17f1c8e584b
ŷ(x1, x2) = sign([x1, x2]' * w + b)

# ╔═╡ 16448a5a-090f-4a6a-a9a0-55192a77dfb3
heatmap!(p, 5:0.1:20, 5:0.1:20, ŷ, color=:coolwarm)

# ╔═╡ Cell order:
# ╟─28042615-8f2b-41fa-9b14-28d69180debb
# ╠═3ec6aae0-c233-11eb-355d-3dfe780e33bc
# ╠═a6997248-aee3-4ebb-8a1e-f791ad503aa8
# ╟─bf7d0459-4b47-4839-800d-01bdcd500633
# ╟─27dfce4e-3ccd-4221-8eca-144ae87e2998
# ╠═5056cc0a-f832-47b3-a6a5-f7ee2d7bc5ba
# ╠═e5334dab-56b2-4dd0-a12c-e17f1c8e584b
# ╠═16448a5a-090f-4a6a-a9a0-55192a77dfb3


text:
@testset "Display Tests" begin
    @testset "AnovaData" begin
        data = [AnovaValue( "Total", 1827.6975, 19),
                AnovaResult(    "A",   70.3125,  1,   70.3125,  3.0706495,  0.098856175),
                AnovaResult(    "B", 1386.1125,  1, 1386.1125, 60.533556,   7.9430782e-7),
                AnovaResult("A × B",    4.9005,  1,    4.9005,  0.21401199, 0.64987001),
                AnovaFactor(    "C",  366.372,  16,   22.89825)]
        result = AnovaData(data, AnovaFactor[], AnovaFactor[], Float64[], 1)
        expectedlines = ["",
                         "Analysis of Variance Results",
                         "",
                         "Effect         SS  DF         MS          F           p",
                         "-------------------------------------------------------",
                         " Total  1827.7     19                                  ",
                         "     A    70.3125   1    70.3125   3.07065   0.0988562 ",
                         "     B  1386.11     1  1386.11    60.5336    7.94308e-7",
                         " A × B     4.9005   1     4.9005   0.214012  0.64987   ",
                         "     C   366.372   16    22.8983                       ",
                         ""]

        @test sprint(show, result) == join(expectedlines, "\n")
    end
end


text:
using BenchmarkTools
using ABM

@btime sin(1)


text:
module Tensors

"A multi-dimensional array decorated with index names"
type Tensor{T,N,Indices<:Tuple} <: AbstractArray{T,N}
    data::Array{T,N}
end

export Tensor

# Conversion from Array should perform some basic checks
import Base.convert
@generated function convert{T,N,Indices<:Tuple}(::Type{Tensor{T,N,Indices}}, data::Array{T,N})
    # Test if Indices has correct length
	if (N == length(Indices.parameters))
		# Test if there are repeated indices
		if (N == length(unique(Indices.parameters)))
			return :(Tensor{$T,$N,$Indices}(data))
		else # We could perform a partial trace, but convert() should preserve data as best as possible.
		    return :(error("Indices $(Indices) are not unique"))
		end
	else
		return :(error("Attempted to assign $(length(Indices.parameters)) indices for $N-dimensional tensor"))
	end
end

# Simple conversion back to an array
convert{T,N,Indices<:Tuple}(::Type{Array{T,N}}, tensor::Tensor{T,N,Indices}) = tensor.data

# Inherit Array definitions and behaviours
import Base.getindex
import Base.setindex!
import Base.eltype
import Base.length
import Base.ndims
import Base.size
import Base.eachindex
import Base.stride
import Base.strides
import Base.linearindexing
import Base.*
import Base.+
import Base.-
import Base.conj
import Base.(.*)
import Base.(./)
import Base.(.^)
import Base.exp
import Base.log
import Base.sin
import Base.cos
import Base.tan
import Base.sinh
import Base.cosh
import Base.tanh
import Base.eig
import Base.svd
import Base.eigs
import Base.svds

# Informational...
eltype{T}(::Tensor{T}) = T
length(t::Tensor) = length(t.data)
ndims{T,N}(::Tensor{T,N}) = N
size(t::Tensor) = size(t.data)
size(t::Tensor,n::Integer) = size(t.data,n)
eachindex(t::Tensor) = eachindex(t.data)
stride(t::Tensor,k::Integer) = eachindex(t.data,k)
strides(t::Tensor) = eachindex(t.data)

# Requirements of AbstractArray (could do DenseArray, I suppose?)
linearindexing{T,N,Indices}(::Type{Tensor{T,N,Indices}}) = linearindexing(Array) # == Base.linearfast()
getindex{T,N,Indices}(t::Tensor{T,N,Indices},i::Int) = getindex(t.data,i)
setindex!{T,N,Indices}(t::Tensor{T,N,Indices},v,i::Int) = setindex!(t.data,v,i)
# similar()...

# Constructing a Tensor from an Array.
@generated function getindex{T,N,Indices<:Tuple}(data::Array{T,N},::Type{Indices})
    idx = svec_to_vec(Indices.parameters)

    # Test if Indices has correct length
	if (N != length(idx)); return :(error("Attempted to assign $(length(Indices.parameters)) indices for $N-dimensional tensor")); end

	# Test if there are repeated indices
	if (N == length(unique(idx))); return :(Tensor{$T,$N,$Indices}(data)); end

	# Perform a (full or partial) trace
	uniques = unique(idx)
    ex_str = "begin; data_tmp = zeros($T,("
	n_sums = 0
	unique_idx = Any[]
	idx_strs = Vector{String}(N)
	for i = 1:length(uniques)
	    x = find(x->x==uniques[i],idx)
	    if length(x) == 1
	        idx_strs[x] = ":"
	        push!(unique_idx,uniques[i])
	        if length(unique_idx) > 1; ex_str = ex_str * ","; end
	        ex_str = ex_str * "size(data,$(x[1]))"
	    end
	end
	ex_str = ex_str * ")); "

	if length(unique_idx) == 0; ex_str = "begin; data_tmp = zero($T); "; end

	for i = 1:length(uniques)
	    x = find(x->x==uniques[i],idx)
	    if length(x) > 1
	        n_sums += 1
	        idx_strs[x] = "i$(n_sums)"
	        ex_str = ex_str * "for i$(n_sums) = 1:length(size(data,$(x[1]))); "
	    end
	end

    ex_str = ex_str * "data_tmp"
    if length(unique_idx) > 0; ex_str = ex_str * "[:]"; end
    ex_str = ex_str * " += data["
    for i = 1:N
        if i > 1; ex_str = ex_str * ","; end
        ex_str = ex_str * "$(idx_strs[i])"
    end
    ex_str = ex_str * "]"
    if length(unique_idx) > 0; ex_str = ex_str * "[:]"; end
    ex_str = ex_str * "; "

	for i = 1:n_sums
	    ex_str = ex_str * "end; "
	end

	if length(unique_idx) > 0
		ex_str = ex_str * "Tensor{$T,$(length(unique_idx)),Tuple{"
		for i = 1:length(unique_idx)
			if i > 1; ex_str = ex_str * ","; end
			ex_str = ex_str * "$(unique_idx[i])"
		end
		ex_str = ex_str * "}}(data_tmp); end"
	else
	    ex_str = ex_str * "data_tmp; end"
	end

	return parse("$ex_str")
end

# Return an Array from a Tensor, permuting indices to put in request ordering
@generated function getindex{T,N,Indices<:Tuple,Out<:Tuple}(tensor::Tensor{T,N,Indices},::Type{Out})
    if  Out == Indices # no permutation necessary
        return :(tensor.data)
    else # Find the permutation
        idx = Indices.parameters
        out = Out.parameters
        if length(out) != N; return :(error("Attempting to de-index a rank-$N tensor with $(length(out)) indices")); end
        if length(out) != length(unique(out)); return :(error("Indices $Out are not unique")); end
        if length(idx) != N || length(idx) != length(unique(idx)); return :(error("Rank-$N tensor has malformed indices $Indices")); end

        return :(permutedims(tensor.data,$(permutation(idx,out))))
    end
end

# Sets an Array given a tensor and order
function setindex!{T,N,Indices<:Tuple,Out<:Tuple}(arrayout::Array{T,N},tensor::Tensor{T,N,Indices},::Type{Out})
    arrayout = tensor[Out]
end

# Multiplication by a scalar
*{Tx<:Number,T,N,Indices}(x::Tx,tensor::Tensor{T,N,Indices}) = Tensor{promote_type(Tx,T),N,Indices}(x*tensor.data)
*{Tx<:Number,T,N,Indices}(tensor::Tensor{T,N,Indices},x::Tx) = Tensor{promote_type(T,Tx),N,Indices}(tensor.data*x)

# Element-wise addition
@generated function +{T1,N,Indices1,T2,Indices2}(tensor1::Tensor{T1,N,Indices1},tensor2::Tensor{T2,N,Indices2})
    if Indices1 == Indices2
        return :(Tensor{$(promote_type(T1,T2)),N,Indices1}(tensor1.data + tensor2.data))
    else # permute tensor2 first
        idx1 = svec_to_vec(Indices1.parameters)
        idx2 = svec_to_vec(Indices2.parameters)
        return :(Tensor{$(promote_type(T1,T2)),N,Indices1}(tensor1.data + permutedims(tensor2.data,$(permutation(idx1,idx2)))))
    end
end

# Subtraction
-{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(-tensor.data)
@generated function -{T1,N,Indices1,T2,Indices2}(tensor1::Tensor{T1,N,Indices1},tensor2::Tensor{T2,N,Indices2})
    if Indices1 == Indices2
        return :(Tensor{$(promote_type(T1,T2)),N,Indices1}(tensor1.data - tensor2.data))
    else # permute tensor2 first
        idx1 = svec_to_vec(Indices1.parameters)
        idx2 = svec_to_vec(Indices2.parameters)
        return :(Tensor{$(promote_type(T1,T2)),N,Indices1}(tensor1.data - permutedims(tensor2.data,$(permutation(idx1,idx2)))))
    end
end

# Element-wise multiplication
@generated function .*{T1,N,Indices1,T2,Indices2}(tensor1::Tensor{T1,N,Indices1},tensor2::Tensor{T2,N,Indices2})
    if Indices1 == Indices2
        return :(Tensor{$(promote_type(T1,T2)),N,Indices1}(tensor1.data .* tensor2.data))
    else # permute tensor2 first
        idx1 = svec_to_vec(Indices1.parameters)
        idx2 = svec_to_vec(Indices2.parameters)
        return :(Tensor{$(promote_type(T1,T2)),N,Indices1}(tensor1.data .* permutedims(tensor2.data,$(permutation(idx1,idx2)))))
    end
end

# Element-wise division
@generated function ./{T1,N,Indices1,T2,Indices2}(tensor1::Tensor{T1,N,Indices1},tensor2::Tensor{T2,N,Indices2})
    if Indices1 == Indices2
        return :(Tensor{$(promote_type(T1,T2)),N,Indices1}(tensor1.data ./ tensor2.data))
    else # permute tensor2 first
        idx1 = svec_to_vec(Indices1.parameters)
        idx2 = svec_to_vec(Indices2.parameters)
        return :(Tensor{$(promote_type(T1,T2)),N,Indices1}(tensor1.data ./ permutedims(tensor2.data,$(permutation(idx1,idx2)))))
    end
end

# Element-wise power
.^{Tx<:Number,T,N,Indices}(tensor::Tensor{T,N,Indices},x::Tx) = Tensor{promote_type(T,Tx),N,Indices}(tensor.data .^ x)

# Common element-wise functions such as exp, log, sin, cos, etc
exp{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(exp(tensor.data))
log{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(log(tensor.data))
sin{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(sin(tensor.data))
cos{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(cos(tensor.data))
tan{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(tan(tensor.data))
sinh{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(sinh(tensor.data))
cosh{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(cosh(tensor.data))
tanh{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(tanh(tensor.data))
sqrt{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(sqrt(tensor.data))

# Complex conjugation
conj{T<:Real,N,Indices}(tensor::Tensor{T,N,Indices}) = tensor
conj{T,N,Indices}(tensor::Tensor{T,N,Indices}) = Tensor{T,N,Indices}(conj(tensor.data))

# Tensor multiply that contracts matching indices
@generated function *{T1,N1,Indices1,T2,N2,Indices2}(tensor1::Tensor{T1,N1,Indices1},tensor2::Tensor{T2,N2,Indices2})
    ex_str = ""

    T = promote_type(T1,T2)
    indices_a = svec_to_vec(Indices1.parameters)
    indices_b = svec_to_vec(Indices2.parameters)

    if length(indices_a) != N1 || length(indices_a) != length(unique(indices_a)); return :(error("Rank-$N tensor has malformed indices $indices_a")); end
    if length(indices_b) != N2 || length(indices_b) != length(unique(indices_b)); return :(error("Rank-$N tensor has malformed indices $indices_b")); end

    ## Prepare the permute order, etc
    #N1 = length(idx1)
    #N2 = length(idx2)
    #Nc = length(intersect(idx1,idx2))
    #(permute1, permute2) = joint_permutation(idx1,idx2)

    # If p

    # Prepare the permute order, etc
    Na = length(indices_a)
    Nb = length(indices_b)
    Nc = length(intersect(indices_a,indices_b))
    Naa = Na - Nc
    Nbb = Nb - Nc
    Nout = Na + Nb - 2*Nc

    permute_aa = zeros(Int64,Na - Nc)
    permute_ac = zeros(Int64,Nc)
    permute_bb = zeros(Int64,Nb - Nc)
    permute_bc = zeros(Int64,Nc)

    iaa = 1
    ibb = 1
    iac = 1
    for i = 1:Na
        # TODO: detect more than two indices... (only occurs if user intervenes)
        if in(indices_a[i],indices_b)
            permute_ac[iac] = i
            iac += 1
        else
            permute_aa[iaa] = i
            iaa += 1
        end
    end
    for i = 1:Nb
        if in(indices_b[i],indices_a)
            # Need it to match perfectly to the entry in permute_ac
            entry_a = find(x->x==indices_b[i],indices_a)[1]
            entry_ac = find(x->x==entry_a,permute_ac)[1]

            permute_bc[entry_ac] = i
        else
            permute_bb[ibb] = i
            ibb += 1
        end
    end

    # Build a string of the multiplication code
    sizes_out = "("
    for i = 1:Naa
        if i != 1
            sizes_out = sizes_out * ","
        end
        sizes_out = sizes_out * "size(tensor1.data,$(permute_aa[i]))"
    end
    for i = 1:Nbb
        if i != 1 || Naa > 0
            sizes_out = sizes_out * ","
        end
        sizes_out = sizes_out * "size(tensor2.data,$(permute_bb[i]))"
    end
    sizes_out = sizes_out * ")"

    if Nout == 0
        # Inner product
        ex_str = "return (tensor1.data[:].' * permutedims(tensor2.data,$(permute_bc))[:])[1]"
    elseif Nc == 0
        # Outer product
        ex_str = "return Tensor{$T,$Nout,Tuple{"
		for i = 1:Naa
		    if i > 1; ex_str = ex_str * ","; end
		    ex_str = ex_str * "$(indices_a[permute_aa[i]])"
		end
		for i = 1:Nbb
		    if i > 1 || Naa > 0; ex_str = ex_str * ","; end
		    ex_str = ex_str * "$(indices_b[permute_bb[i]])"
		end
		ex_str = ex_str * "}}(reshape(tensor1.data[:] * tensor2.data[:].',$sizes_out))"
    else
        # Matrix product (or matrix-vector/vector-matrix)

        size_aa = "1"
        for i = 1:Naa
            size_aa = size_aa * "*size(tensor2.data,$(permute_aa[i]))"
        end

        size_c = "1"
        for i = 1:Nc
            size_c = size_c * "*size(tensor1.data,$(permute_ac[i]))"
        end

        size_bb = "1"
        for i = 1:Nbb
            size_bb = size_bb * "*size(tensor2.data,$(permute_bb[i]))"
        end

		ex_str = "return Tensor{$T,$Nout,Tuple{"
		for i = 1:Naa
		    if i > 1; ex_str = ex_str * ","; end
		    ex_str = ex_str * "$(indices_a[permute_aa[i]])"
		end
		for i = 1:Nbb
		    if i > 1 || Naa > 0; ex_str = ex_str * ","; end
		    ex_str = ex_str * "$(indices_b[permute_bb[i]])"
		end
		ex_str = ex_str * "}}(reshape(reshape(permutedims(tensor1.data,$(vcat(permute_aa,permute_ac))),($size_aa,$size_c)) * reshape(permutedims(tensor2.data,$(vcat(permute_bc,permute_bb))),($size_c,$size_bb)),$sizes_out))"
    end

    return parse(ex_str)
end

# Eigenvalue with indices arranged
@generated function eig{T,N,Indices,IndicesLeft<:Tuple,IndicesRight<:Tuple,IndicesEig<:Tuple}(tensor::Tensor{T,N,Indices},::Type{IndicesLeft},::Type{IndicesRight},index_eig::Type{IndicesEig} = Tuple{:eig})
    # We need to permute the matrix into order [indices_left, indices_right]
    idx = svec_to_vec(Indices.parameters)
    idx_left = svec_to_vec(IndicesLeft.parameters)
    idx_right = svec_to_vec(IndicesRight.parameters)
    idx_permuted = vcat(idx_left,idx_right)
    idx_eig = svec_to_vec(IndicesEig.parameters)

    if length(idx) != length(idx_permuted); return :(error("Number of indices ($(length(idx_left)), $(length(idx_right))) does not match tensor's $N")); end
    if length(idx_left) == 0; return :(error("At least one left index must be specified")); end
    if length(idx_right) == 0; return :(error("At least one right index must be specified")); end
    if length(idx_eig) != 1; return :(error("Only one index should be associated with the eigenvalues")); end

    p = permutation(idx,idx_permuted)

    # sizes
    size_left = "size(tensor.data,$(p[1]))"
    resize_left = "size(tensor.data,$(p[1]))"
    for i = 2:length(idx_left)
        size_left = size_left * "*size(tensor.data,$(p[i]))"
        resize_left = resize_left * ",size(tensor.data,$(p[i]))"
    end
    size_right = "size(tensor.data,$(p[1+length(idx_left)]))"
    #resize_right = "size(tensor.data,$(p[1+length(idx_left)]))"
    for i = 2:length(idx_right)
        size_right = size_right * "*size(tensor.data,$(p[i+length(idx_left)]))"
    #    resize_right = resize_right * ",size(tensor.data,$(p[i+length(idx_left)]))"
    end

    str_left = Vector{String}(length(idx_left))
    for i = 1:length(idx_left); isa(idx_left[i],Symbol) ? str_left[i] = ":$(idx_left[i])" : str_left[i] = "$(idx_left[i])"; end
    #str_right = Vector{String}(length(idx_right))
    #for i = 1:length(idx_right); isa(idx_right[i],Symbol) ? str_right[i] = ":$(idx_right[i])" : str_left[i] = "$(idx_right[i])"; end
    str_eig = (isa(idx_eig[1],Symbol) ? ":$(idx_eig[1])" : "$(idx_eig[1])")

    ex_str = "begin; "
    if idx_permuted == idx # No permutation necessary, only reshape.
        ex_str = ex_str * "(D,V) = eig(reshape(tensor.data,($size_left,$size_right))); "
    else # permutation is necessary
        ex_str = ex_str * "(D,V) = eig(reshape(permutedims(tensor.data,$p),($size_left,$size_right))); "
    end

	ex_str = ex_str * "(Tensor{$T,1,Tuple{$(str_eig)}}(D),Tensor{$T,$(1+length(idx_left)),Tuple{"
	for i = 1:length(idx_left)
		ex_str = ex_str * "$(str_left[i]),"
	end
	ex_str = ex_str * "$(str_eig)}}(reshape(V,($(resize_left),$(size_left))))); "

    ex_str = ex_str * "end"
    #display(ex_str)
    return parse(ex_str)
end

@generated function svd{T,N,Indices,IndicesLeft<:Tuple,IndicesRight<:Tuple,IndicesSVD<:Tuple}(tensor::Tensor{T,N,Indices},::Type{IndicesLeft},::Type{IndicesRight},index_svd::Type{IndicesSVD} = Tuple{:svd})
    # We need to permute the matrix into order [indices_left, indices_right]
    idx = svec_to_vec(Indices.parameters)
    idx_left = svec_to_vec(IndicesLeft.parameters)
    idx_right = svec_to_vec(IndicesRight.parameters)
    idx_permuted = vcat(idx_left,idx_right)
    idx_svd = svec_to_vec(IndicesSVD.parameters)

    if length(idx) != length(idx_permuted); return :(error("Number of indices ($(length(idx_left)), $(length(idx_right))) does not match tensor's $N")); end
    if length(idx_left) == 0; return :(begin; n = vecnorm(tensor.data); (one(T),n,tensor/n); end); end
    if length(idx_right) == 0; return :(begin; n = vecnorm(tensor.data); (tensor/n,n,one(T)); end); end
    if length(idx_svd) != 1; return :(error("Only one index should be associated with the singular values")); end

    p = permutation(idx,idx_permuted)

    # sizes
    size_left = "size(tensor.data,$(p[1]))"
    resize_left = "size(tensor.data,$(p[1]))"
    for i = 2:length(idx_left)
        size_left = size_left * "*size(tensor.data,$(p[i]))"
        resize_left = resize_left * ",size(tensor.data,$(p[i]))"
    end
    size_right = "size(tensor.data,$(p[1+length(idx_left)]))"
    resize_right = "size(tensor.data,$(p[1+length(idx_left)]))"
    for i = 2:length(idx_right)
        size_right = size_right * "*size(tensor.data,$(p[i+length(idx_left)]))"
        resize_right = resize_right * ",size(tensor.data,$(p[i+length(idx_left)]))"
    end

    str_left = Vector{String}(length(idx_left))
    for i = 1:length(idx_left); isa(idx_left[i],Symbol) ? str_left[i] = ":$(idx_left[i])" : str_left[i] = "$(idx_left[i])"; end
    str_right = Vector{String}(length(idx_right))
    for i = 1:length(idx_right); isa(idx_right[i],Symbol) ? str_right[i] = ":$(idx_right[i])" : str_right[i] = "$(idx_right[i])"; end
    str_svd = (isa(idx_svd[1],Symbol) ? ":$(idx_svd[1])" : "$(idx_svd[1])")

    ex_str = "begin; "
    if idx_permuted == idx # No permutation necessary, only reshape.
        ex_str = ex_str * "(U,S,V) = svd(reshape(tensor.data,($size_left,$size_right))); "
    else # permutation is necessary
        ex_str = ex_str * "(U,S,V) = svd(reshape(permutedims(tensor.data,$p),($size_left,$size_right))); "
    end

    ex_str = ex_str * "(Tensor{$T,$(1+length(idx_left)),Tuple{"
	for i = 1:length(idx_left)
		ex_str = ex_str * "$(str_left[i]),"
	end
	ex_str = ex_str * "$(str_svd)}}(reshape(U,($(resize_left),$(size_left)))), Tensor{$T,1,Tuple{$(str_svd)}}(S), Tensor{$T,$(1+length(idx_right)),Tuple{"
	for i = 1:length(idx_right)
		ex_str = ex_str * "$(str_right[i]),"
	end
	ex_str = ex_str * "$(str_svd)}}(reshape(V,($(resize_right),$(size_right))))); "

    ex_str = ex_str * "end"
    display(ex_str)
    return parse(ex_str)
end

#################
### Utilities ###
#################

"Convert Simple Vector to Vector"
function svec_to_vec(svec::SimpleVector)
    vec = Vector{Any}(length(svec))
    for i = 1:length(svec)
        vec[i] = svec[i]
    end
    return vec
end

"Function finds the permutation that transforms in into out"
function permutation(in,out)
    N = length(in)
    if length(out) != N; error("Permutation must be between equal-sized sets"); end
    permutation = zeros(Int64,N)
    for i = 1:N
        m = 0
        matches = 0
        for j = 1:N
            if in[i] == out[j]
                matches += 1
                m = j
            end
        end
        if matches != 1
            error("Output $out must match up with input $in")
        end
        permutation[i] = m
    end
    return permutation
end

##############
### Macros ###
##############
"Macro takes expression and replaces [[ a, ... ]] with [Tuple{:a,...}]"
macro tensor(ex)
    # search through ex and replace all instances of vect with something else
    replace_ref_vect_with_tuple!(ex)

    return esc(ex) # esc() means ex gets returned as-is for running in the caller's scope
end

export @tensor

#######################
### Code generation ###
#######################
"Recursive function that replaces [[ a, ... ]] with [Tuple{:a,...}] in an expression tree"
function replace_ref_vect_with_tuple!(a::Expr)
    for i = 1:length(a.args)
		if isa(a.args[i],Expr)
		    replace_ref_vect_with_tuple!(a.args[i])
		end
	end

    if a.head == :ref
        if isa(a.args[2],Expr) && a.args[2].head == :vect
            # BINGO - we found the double bracket [[ ... ]]

            tmp = Vector{Any}(1+length(a.args[2].args))
            tmp[1] = :Tuple
            for i = 1:length(a.args[2].args)
                tmp2 = (a.args[2].args[i])
                if isa(tmp2,Symbol)
                    tmp[i+1] = :($tmp2)
                else
                    tmp[i+1] = tmp2
                end
            end
            a.args[2].head = :curly
            a.args[2].args = tmp
        end
    end
end

end # module


text:
module Clusterpath

using LinearAlgebra
using Random
using DataFrames
using Interpolations
using Distributed
using MultivariateStats
using Distributions
import Gaston

include("AVLTree.jl")

export clusterpath, plot_cluster, plot_path, assign_cluster
export clusterpath_pop, cond_mean_on_LR, find_truncation, find_split
export generate_mixture_normal

"""
    clusterpath!(y; α=0.1, jitter=1e-6, return_split=false, silence_warning=false)
Perform convex clustering via α-thresholded clusterpath algorithm to `y`.
To avoid ties, add a small perturbation (`jitter`) to each data, and then sort it.
Returns a dictionary of
* (1) idx: cluster indices assigned to each elements in `y`
* (2) lambda: λ values in each merge
* (3) alpha: solution corresponding to each observation where merge happened
* (4) n_cluster: number of clusters wrt merging algorithm
* (5) n_splits: number of total splits, in corresponding top-down algorithm
* (6) splits: sample split point (only get returned when `return_split=true`)
"""
function clusterpath(y; α=0.1, jitter=1e-6, return_split=false, silence_warning=false)
    n = length(y)

    # remove ties with small perturbation `jitter`.
    x = copy(y)
    Random.seed!(0)
    x += randn(n) * jitter

    ii = sortperm(sortperm(x))
    sort!(x)  # sort `x` in ascending order

    # Initialize variables
    n_cluster = n  # (n) number of clusters
    probability = 0

    ## linkage between clusters
    next = Array{Union{Missing,Int64},1}(2:n+1)
    next[n] = missing
    prev = Array{Union{Missing,Int64},1}(0:n-1)
    prev[1] = missing

    ## status of clusters
    cardinals = ones(Int64, n)  # (freq) cardinality(freq) of each clusters
    centroids = copy(x)
    distance = diff(centroids) ./ (cardinals[1:n_cluster-1] + cardinals[2:end])
    λ_old = 0
    idx_cluster = Array(1:n) # (ind) cluster indices of observations

    ## big merge tracker
    was_big = false  # true if last merge was big
    is_big = false   # true if current merge is big

    # store returns along the solution path
    alphas = Array{Array{Float64,1},1}([copy(centroids)[ii]])
    λs = Float64[λ_old]  # (lambda) lambdas
    idx_along_path = Array{Array{Int64,1},1}()#[Array(1:n)[ii]])
    n_cluster_along_path = Int[]#n]
    n_splits = Int64(0)
    splits = Float64[]
#    println(idx_cluster[ii])
#    println(round.(distance,digits=5))

    # tree initailization
    ## (tree of pairwise distance btw successive indices)
    tree = AVLTree{Float64}()
    for i in 1:n-1
        push!(tree, distance[i], i)
    end

    """
    function for status update at each lambda, except for alpha
    """
    function update_return!(; update_lambda=true)
        if update_lambda
            append!(λs, λ_old)
        end
        append!(idx_along_path, [idx_cluster[ii]])
        append!(n_cluster_along_path, n_cluster)
    end

    """
    function for status update just for alpha
    """
    function update_alpha!()
        walk = 1
        alpha = Array{Float64}(undef, n)
        while !ismissing(walk)
            if !ismissing(next[walk])
                #n_repeat = next[walk] - walk
                alpha[walk:next[walk]-1] .= centroids[walk] + λ_old *((n-next[walk]+1)-(walk-1))
            else
                alpha[walk:end] .= centroids[walk] + λ_old *(-(walk-1))
            end
            walk = next[walk]
        end
        append!(alphas, [alpha[ii]])
    end

    ## main loop
    while n_cluster > 1

        # find minimum distance
        pivot = minimum_node(tree.root)
        λ_old, i = pivot.data, pivot.idx  # λ of next merge (minimal distance)

        #println(cardinals[i], " ", cardinals[next[i]])
#        println((i, next[i]))

        # detect and handle big merge
        if min(cardinals[i], cardinals[next[i]]) > ceil(α * n)  # But the paper says ⌈n * α⌉
            #println(cardinals[i], " ", cardinals[next[i]], " ", i, " ", λ_old, " ", n)

            is_big = true
            probability = (cardinals[i] + cardinals[next[i]]) / n

            n_splits += 1
            if return_split
                a_1 = max((x[idx_cluster .== idx_cluster[i]])...)
                a_2 = min((x[idx_cluster .== idx_cluster[next[i]]])...)
                split = (cardinals[next[i]] * a_2 + cardinals[i] * a_1) / (cardinals[next[i]] + cardinals[i])
                append!(splits, split)
            end
        end

        if is_big || was_big
            if is_big
                was_big = true
                is_big = false
            else
                was_big = false
            end
            # update λ path
            update_return!()
        end

        # update status
        centroids[i] = (centroids[i] * cardinals[i] + centroids[next[i]] * cardinals[next[i]]) / (cardinals[i] + cardinals[next[i]])
        cardinals[i] += cardinals[next[i]]
        idx_cluster[idx_cluster .== next[i]] .= i
        #n_right[i] = n_right[next[i]]

        # update distance and tree
        delete!(tree, distance[i])

        if !ismissing(prev[i])
            delete!(tree, distance[prev[i]])
            d_prev = (centroids[i] - centroids[prev[i]]) / (cardinals[prev[i]] + cardinals[i])
            distance[prev[i]] = d_prev
            push!(tree, d_prev, prev[i])
        end
        if !ismissing(next[next[i]])
            delete!(tree, distance[next[i]])
            d_next = (centroids[next[next[i]]] - centroids[i]) / (cardinals[i] + cardinals[next[next[i]]])
            distance[i] = d_next
            push!(tree, d_next, i)
        end

        # update linkage
        next[i] = next[next[i]]
        if !ismissing(next[i])
            prev[next[i]] = i
        end

        n_cluster -= 1

        if is_big || was_big
            if is_big
                was_big = true
                is_big = false
            else
                was_big = false
            end
            # update λ path
            #update_return!()
            update_alpha!()
        end
    end

    # last update for status at n_cluster = 1
    update_return!(update_lambda=false)

    # if the last Big Merge has sum of probabilities < 50%,
    # consider the case as if no splits have ever occured
    if probability < 0.5
        if !silence_warning
            @warn "No clusters were found: the last big merge has probability $(probability) < 0.5"
        end
        idx_along_path = Array(1:n)
        n_cluster_along_path = n
        n_splits = 0
    end

    return Dict("lambda" => λs,
                "idx" => idx_along_path,
                "alpha" => alphas,
                "n_cluster" => n_cluster_along_path,
                "n_splits" => n_splits,
                "splits" => splits)

end

"""
    rep(x, times)
Julia equivalent of R's `rep`.
"""
function rep(x, times)
    vcat(fill.(x, times)...)
end

"""
    merge_dfs(df1, df2)
merge two solution path dataframes with `outerjoin()`
"""
function merge_dfs(df1, df2)
    outerjoin(df1, df2; on=[:obs, :lambda], makeunique=true)
end

"""
    cast_solution(x; α=0.1, silence_warning=false)

A generator function that casts a dataframe of solution path using `clusterpath()`. The function yields the solution path of each observation one at an iteration. This return will be used in `plot_path()` and `plot_cluster()`.

* `x`: the data for which we want to cluster.
* `α`: threshold for Big Merge Tracker procedure. `α` shoulbe be ≥ 0. The smaller the value of `α`, the more detailed the solution path will be. The larger the `α`, the faster the algorithm will run. (default: `0.1`)
"""
function cast_solution(x; α=0.1, silence_warning=false)

    n, p = try
             n, p = size(x)
             (n, p)
           catch e
             if isa(e, BoundsError)
               (length(x), 1)
             else
               throw(e)
             end
           end
#     println(n, " ", p)

    #lambdas = Array{Union{Nothing,Array{Float64,1}},1}(repeat([nothing], p))
    ccs = Array{Dict,1}(undef, p)
    lambdas = Array{Array{Float64,1},1}(undef, p)
    Threads.@threads for k=1:p
        cc = clusterpath(x[:, k]; α=α, silence_warning=silence_warning)
        ccs[k] = cc
        lambdas[k] = cc["lambda"]
    end

    lambdas = vcat(lambdas...)
    unique!(lambdas)
    sort!(lambdas)

    Channel{DataFrame}() do channel
        # report one obs at a time (for memory management)
        for i=1:n
            # initialize solution path matrix for ith obs
            solution = Array{Union{Float64, Missing}, 2}(missing, length(lambdas), p+2)
            solution[:, 1] .= i
            solution[:, 2] = lambdas

            # loop over dimensions
            for k=1:p
                cc = ccs[k]

                idx = [l in cc["lambda"] for l in lambdas]

                alpha_idx = 1
                for (lam_idx, is_in_lambda) in enumerate(idx)
                    if is_in_lambda
                        solution[lam_idx, k+2] = cc["alpha"][alpha_idx][i]
                        alpha_idx += 1
                    end
                end

#                 alpha = Float64[a[i] for a in cc["alpha"]]

#                 df = DataFrame(obs = repeat([i], n_path),
#                                lambda = lambdas[k],
#                                dim = alpha)#,
#                                #col = repeat(["dim$k"], n * n_path))

#                 if isnothing(solution)
#                     solution = df
#                 else
#                     solution = merge_dfs(solution, df)
#                 end

                # calculate cluster number based on n_splits
            end

            solution = DataFrame(solution)
            rename!(solution, vcat(["obs", "lambda"], ["dim$k" for k=1:p]))
#             sort!(solution, [:obs, :lambda])

            put!(channel, solution)
        end
    end
end

"""
    _prepare(solution_i; for_path::Bool)

prepare the return of `cast_solution()` for `plot_path()`, or `plot_cluster()`.
* `solution_i`: single solution path element of `cast_solution()`.
* `for_path`: if `true` then prepare for `plot_path()`, `false` then prepare for `plot_cluster()`.
"""
function _prepare!(solution_i; for_path::Bool)
    p = size(solution_i, 2) - 2
    n = length(unique(solution_i.obs))

    #unique!(solution_i, vcat(["obs"], ["dim$k" for k=1:p]))
    #Impute.interp!(solution_i)

    if !for_path
        nothing
        #solution_i[1:end-n_node-1, :] .= missing
        #dropmissing!(solution_i)
    end

    for j=1:p
        idx_missing = ismissing.(solution_i["dim$j"])
        #trunc = findlast(.~idx_missing)

        lambda_obs = solution_i[.~idx_missing, "lambda"]
        point_obs = Float64.(solution_i[.~idx_missing, "dim$j"])
        itp = interpolate((lambda_obs,), point_obs, Gridded(Linear()))
        etp = extrapolate(itp, Flat())

        solution_i[idx_missing, "dim$j"] = etp.(solution_i.lambda[idx_missing])
    end

    #Impute.locf!(solution_i)
end

"""
    plot_path(x; α=0.1, gt=ones(Int64, size(x,1)), force2dim=true, show::Bool, savefig=false, fname="path_plot", verbose=false)

plot clusterpath with the data(`x`).
If the dimension of `x` is greater than 2, only plot a combination of first two dimensions.

***`Gaston.jl` and `gnuplot` should be installed and on the PATH of your system. Install gnuplot [here](https://sourceforge.net/projects/gnuplot/files/gnuplot/).***

* `x`: the data to cluster.
* `α`: threshold for BMT-clusterpath
* `gt`: the ground truth labels.
* `force2dim`: force dim>2 data to 2-dim data (default: `true`)
* `show`: whether to show the plot in the notebook. **Highly recommended not to show if the number of samples is large.**
* `savefig`: whether to save the figure as a PNG file. (default: `false`)
* `fname`: image file name to be used when `savefig` is `true`. (default: `"path_plot"`)
"""
function plot_path(x; α=0.1, gt=ones(Int64, size(x,1)), force2dim=true, show::Bool, savefig=false, fname="path_plot", verbose=false)
    n_dim = length(size(x)) == 1 ? 1 : size(x, 2)
    if n_dim > 2 && !force2dim
        x = x[:, 1:3]
    elseif (n_dim > 2 && force2dim) || n_dim == 2
        x = x[:, 1:2]
    end
    solution = cast_solution(x, α=α, silence_warning=true)
#     n_plot = Int64(n_dim * (n_dim-1) / 2)
#     if n_plot > 6
#         n_plot = 6
#         n_dim = 4
#     end

    # initialize plot
#     if n_dim > 2
#         p = Plots.plot(layout=(1,n_plot), size=(400*n_plot, 400), legend=false)
#         counter = 0
#             for i=1:n_dim
#                 for j=(i+1):n_dim
#                     counter += 1
#                     Plots.plot!(p[counter], view(x, :, i), view(x, :, j), seriestype=:scatter, color=gt, alpha=.2)
#                 end
#             end
    if n_dim > 2 && !force2dim
        #p = Plots.plot(view(x, :, 1), view(x, :, 2),
        #               #series_annotations = text.(1:size(x, 1), :bottom),
        #               seriestype=:scatter, color=gt, alpha=.2, legend=false)
        clst = unique(gt)
        p = Gaston.scatter3(
                x[:, 1], x[:, 2], x[:, 3], supp=gt .% 8, variable="using 1:2:3:(4)",
                pointtype="ecircle", lw=0.75, lc="palette",
                Gaston.Axes(key="noautotitle",
                            grid=:on,
                            style="fill transparent solid 0.005;
                                   set style circle radius screen 0.005;
                                   set border 21;
                                   unset colorbox;
                                   set auto fix;
                                   set palette defined ( 0 '#E41A1C', 1 '#377EB8', 2 '#4DAF4A', 3 '#984EA3', 4 '#FF7F00', 5 '#d2d200', 6 '#A65628', 7 '#F781BF' );
                                   set cbrange [0:7];
                                   set size square 1,1;
                                   set format y \"\"; set format x \"\"; set format z \"\"")
            )
#         p = Gaston.plot(
#                 x[:, 1], x[:, 2], z=x[:, 3], supp=gt.%8, variable="using 1:2:3:4", w=:circles, lw=0.75, lc="palette",
#                 Gaston.Axes(key="noautotitle",
#                             style="fill transparent solid 0.005;
#                                    set style circle radius screen 0.005;
#                                    set border 3;
#                                    unset colorbox;
#                                    set auto fix;
#                                    set palette defined ( 0 '#E41A1C', 1 '#377EB8', 2 '#4DAF4A', 3 '#984EA3', 4 '#FF7F00', 5 '#d2d200', 6 '#A65628', 7 '#F781BF' );
#                                    set cbrange [0:7];
#                                    set size square 1,1;
#                                    set format y \"\"; set format x \"\"")
#             )
    elseif (n_dim > 2 && force2dim) || n_dim == 2
        #p = Plots.plot(view(x, :, 1), view(x, :, 2),
        #               #series_annotations = text.(1:size(x, 1), :bottom),
        #               seriestype=:scatter, color=gt, alpha=.2, legend=false)
        clst = unique(gt)
        p = Gaston.plot(
                x[:, 1], x[:, 2], supp=gt.%8, variable="using 1:2:3", w=:circles, lw=0.75, lc="palette",
                Gaston.Axes(key="noautotitle",
                            grid=:on,
                            style="fill transparent solid 0.005;
                                   set style circle radius screen 0.005;
                                   set border 3;
                                   unset colorbox;
                                   set auto fix;
                                   set palette defined ( 0 '#E41A1C', 1 '#377EB8', 2 '#4DAF4A', 3 '#984EA3', 4 '#FF7F00', 5 '#d2d200', 6 '#A65628', 7 '#F781BF' );
                                   set cbrange [0:7];
                                   set size square 1,1;
                                   set format y \"\"; set format x \"\"")
            )
        #if length(clst) > 1
        #    for cl in clst[2:end]
        #        Gaston.plot!(x[gt .== cl, 1], x[gt .== cl, 2], w=:circles, lw=0.75)
        #    end
        #end
    else
        n = size(x, 1)
        #p = Plots.plot(zeros(n), x, seriestype=:scatter, color=gt, alpha=.2, legend=false)

        ## USE GASTON (GNUPLOT) FOR FASTER PLOTTING
        clst = unique(gt)
        p = Gaston.plot(
                zeros(n), x, supp=gt.%8, variable="using 1:2:3", w=:circles, lw=0.75, lc="palette",
                Gaston.Axes(key="noautotitle",
                            grid=:on,
                            style="fill transparent solid 0.005;
                                   set style circle radius screen 0.01;
                                   set border 3;
                                   unset colorbox;
                                   set auto fix;
                                   set palette defined ( 0 '#E41A1C', 1 '#377EB8', 2 '#4DAF4A', 3 '#984EA3', 4 '#FF7F00', 5 '#d2d200', 6 '#A65628', 7 '#F781BF' );
                                   set cbrange [0:7];
                                   set size square 1,1;
                                   set format y \"\"; set format x \"\"")
            )
        #if length(clst) > 1
        #    for cl in clst[2:end]
        #        Gaston.plot!(zeros(n)[gt .== cl], x[gt .== cl], w=:circles, lw=0.75)
        #    end
        #end
    end

    # draw paths, one at a time
    i = 1
    for solution_i in solution
        _prepare!(solution_i; for_path=true)
        #println(solution_i)

#         if n_dim > 2
#             counter = 0
#             for i=1:n_dim
#                 for j=(i+1):n_dim
#                     counter += 1
#                     Plots.plot!(p[counter],
#                                 view(solution_i, "dim$i"), view(solution_i, "dim$j"),
#                                 seriestype=:path, color="black", alpha=0.25, legend=false)
#                 end
#             end
        if n_dim > 2 && !force2dim
            verbose && (i % 500 == 0) && print("$(i) ")

            #Plots.plot!(p, view(solution_i, "dim1"), view(solution_i, "dim2"),
            #            seriestype=:path, color="black", alpha=0.25, legend=false)
            Gaston.scatter3!(view(solution_i, "dim1"), view(solution_i, "dim2"), view(solution_i, "dim3"),
                            w=:lines, lw=0.05, lc="black")
#             Gaston.plot!(view(solution_i, "dim1"), view(solution_i, "dim2"), z=view(solution_i, "dim3"),
#                          w=:lines, lw=0.05, lc="black")
        elseif (n_dim > 2 && force2dim) || n_dim == 2
            verbose && (i % 500 == 0) && print("$(i) ")

            #Plots.plot!(p, view(solution_i, "dim1"), view(solution_i, "dim2"),
            #            seriestype=:path, color="black", alpha=0.25, legend=false)
            Gaston.plot!(view(solution_i, "dim1"), view(solution_i, "dim2"),
                         w=:lines, lw=0.05, lc="black")
        else
            verbose && (i % 500 == 0) && print("$(i) ")

            #Plots.plot!(p, view(solution_i, "lambda"), view(solution_i, "dim1"),
            #            seriestype=:path, color="black", alpha=0.25, legend=false)
            Gaston.plot!(view(solution_i, "lambda"), view(solution_i, "dim1"),
                         w=:lines, lw=0.1, lc="black")
        end

        i += 1
    end

    if savefig
        #png(p, fname*".png")
        Gaston.save(term = "png",
                    output= "$(fname).png",
                    size = "500,500",
                    background = "white")
    end

    if show
        return(p)
    end
end

"""
    assign_cluster(x; α=0.1, verbose=false)

assign cluster to each of the observations in `x`.
returns an array of length=size(x, 1) of cluster indices.

* `x`: data
* `α`: threshold for BMT-clusterpath
"""
function assign_cluster(x; α=0.1, verbose=false)

    n = size(x, 1)
    p = length(size(x)) == 1 ? 1 : size(x, 2)

    final = Array{Float64,2}(undef, n, p)

    clst = repeat([""], n)
    for i=1:p
        splits = clusterpath(x[:, i], α=α, return_split=true, silence_warning=true)["splits"]
        clst_dim = zeros(Int64, n)
        for s in splits
            clst_dim .+= Int64.(x[:, i]  .> s)
        end
        clst .*= string.(clst_dim)
    end

    DataFrames.levelcode.(DataFrames.CategoricalArray(clst))
end
# function assign_cluster!(x, solution; n_node=1, verbose=false)

#     n = size(x, 1)
#     p = length(size(x)) == 1 ? 1 : size(x, 2)

#     final = Array{Float64,2}(undef, n, p)
#     lambda = nothing

#     i = 1
#     for solution_i in solution
#         if verbose && (i % 500 == 0)
#             print("$(i) ")
#         end
#         _prepare!(solution_i; for_path=false, n_node=n_node)

#         #println(solution_i[end, 3:end])
#         final[i, :] = collect(solution_i[end-n_node, 3:end])
#         isnothing(lambda) && (lambda = solution_i.lambda[end-n_node])
#         i += 1
#     end

#     centroids = unique(final, dims=1)
#     n_centroids = size(centroids, 1)

#     idx_x = Array{Int64}(undef, n)

#     Threads.@threads for i=1:n_centroids
#         filter = vec(reduce(&, final .== centroids[i, :]', dims=2))
#         idx_x[filter] .= i
#     end

#     idx_x, lambda
# end

"""
    plot_cluster(x; α=0.1, n_node=1, show::Bool, savefig=false, fname="plot_clst", verbose=false)

Plots the scatter plot of the data `x` colored according to the cluster assigned by clusterpath algorithm.
If the dimension of `x` is greater than 2, perform PCA and plot two PCs.

***`Gaston.jl` and `gnuplot` should be installed and on the PATH of your system. Install gnuplot [here](https://sourceforge.net/projects/gnuplot/files/gnuplot/). ***

* `x`: data
* `α`: threshold for BMT-clusterpath
* `n_node`: if greater than `1`, will assign clusters from previous merge status. (default: `1`)
* `show`: whether to show the figure.
* `savefig`: whether to save the figure as a png file. (default: `false`)
* `fname`: file name to save if `savefig` is true. (default: `"plot_clst"`)
* `verbose`: print out current iteration. (default: `false`)
"""
function plot_cluster(x; α=0.1, force2dim=true, show::Bool, savefig=false, fname="plot_clst", verbose=false)
    clst = assign_cluster(x; α=α, verbose=verbose)

    #solution = cast_solution(x, α=α)
    #solution_1 = take!(solution)

    k = levels(clst)[end]

    p = size(x, 2)

    if p > 3
        force2dim && (pca = fit(PCA, Array(x'); maxoutdim=2))
        !force2dim && (pca = fit(PCA, Array(x'); maxoutdim=3))
        x = MultivariateStats.transform(pca, x')'
    elseif p == 3
        force2dim && (pca = fit(PCA, Array(x'); maxoutdim=2); x = MultivariateStats.transform(pca, x')')
    end

    # use Gaston for plotting
    #uniq_clusters = unique(clst)
    #cl = uniq_clusters[1]
    #idx = clst .== cl
    if p > 2 && !force2dim
        plt = Gaston.scatter3(
                x[:, 1], x[:, 2], x[:, 3], supp=clst .% 8, variable="using 1:2:3:(4)",
                pointtype="ecircle", lw=0.75, lc="palette",
                Gaston.Axes(key="noautotitle",
                            grid=:on,
                            style="fill transparent solid 0.005;
                                   set style circle radius screen 0.005;
                                   set border 21;
                                   unset colorbox;
                                   set auto fix;
                                   set palette defined ( 0 '#E41A1C', 1 '#377EB8', 2 '#4DAF4A', 3 '#984EA3', 4 '#FF7F00', 5 '#d2d200', 6 '#A65628', 7 '#F781BF' );
                                   set cbrange [0:7];
                                   set size square 1,1;
                                   set format y \"\"; set format x \"\"; set format z \"\"")
            )
    elseif (p > 2 && force2dim) || p == 2
        plt = Gaston.plot(
            x[:, 1], x[:, 2], supp=clst .% 8, variable="using 1:2:3",
            w=:circles, lw=0.75, lc="palette",
            Gaston.Axes(key="noautotitle",
                        grid=:on,
                        xlabel="\"k=$(k)\"",
                        style="fill transparent solid 0.005;
                               set style circle radius screen 0.005;
                               set border 3;
                               unset colorbox;
                               set auto fix;
                               set palette defined ( 0 '#E41A1C', 1 '#377EB8', 2 '#4DAF4A', 3 '#984EA3', 4 '#FF7F00', 5 '#d2d200', 6 '#A65628', 7 '#F781BF' );
                               set cbrange [0:7];
                               set size square 1,1;
                               set format y \"\"; set format x \"\"")
        )
    else
        plt = Gaston.plot(
            x, randn(length(x)) * 0.02, supp=clst .% 8, variable="using 1:2:3",
            w=:circles, lw=0.75, lc="palette",
            Gaston.Axes(key="noautotitle",
                        grid=:on,
                        xlabel="\"k=$(k)\"",
                        style="fill transparent solid 0.005;
                               set style circle radius screen 0.005;
                               set border 3;
                               set yrange [-1:1];
                               unset colorbox;
                               set auto fix;
                               set palette defined ( 0 '#E41A1C', 1 '#377EB8', 2 '#4DAF4A', 3 '#984EA3', 4 '#FF7F00', 5 '#d2d200', 6 '#A65628', 7 '#F781BF' );
                               set cbrange [0:7];
                               set size square 1,1;
                               set format y \"\"; set format x \"\"")
        )
    end

#     if length(uniq_clusters) > 1
#         for cl in uniq_clusters[2:end]
#             idx = clst .== cl
# #             println(sum(idx))
#             Gaston.plot!(x[idx, 1], x[idx, 2], w=:circles, lw=0.75)
#         end
#     end

    if savefig
        Gaston.save(term = "png",
                    output= "$(fname).png",
                    size = "500,500",
                    background = "white")
    end

#     plot(x[:, 1], x[:, 2], color=clst,
#          seriestype=:scatter, legend=false, alpha=.2,
#          title="lambda=$(round(λ, digits=3)), k=$(k)")
    if show
        return(plt)
    end
end

function count_final_cluster(x; α=0.1)
    levels(assign_cluster(x, α=α))[end]
end

# ---- population procedure ---- #

"""
    pnorm(x)

Julia equivalent of R's `pnorm`
"""
function pnorm(x)
    cdf(Normal(), x)
end

"""
    dnorm(x)

Julia equivalent of R's `dnorm`
"""
function dnorm(x)
    pdf(Normal(), x)
end

# """
#     cond_mean_on_LR(L::T, R::T, a1::T, a2::T, a3::T, μ1::T, μ2::T, μ3::T) where T <: Real

# Conditional mean on (L, R), defined as μ_{L,R} = (∫_L^R f(x) dx)^(-1) / (∫_L^R x f(x) dx)
# """
# function cond_mean_on_LR(L::T, R::T, a1::T, a2::T, a3::T, μ1::T, μ2::T, μ3::T) where T <: Real
#     pL = a1 * pnorm(L-μ1) + a2 * pnorm(L-μ2) + a3 * pnorm(L-μ3)
#     pR = a1 * pnorm(R-μ1) + a2 * pnorm(R-μ2) + a3 * pnorm(R-μ3)

#     ∫_L = -(a1 * dnorm(L-μ1) + a2 * dnorm(L-μ2) + a3 * dnorm(L-μ3)) +
#           a1 * μ1 * pnorm(L-μ1) + a2 * μ2 * pnorm(L-μ2) + a3 * μ3 * pnorm(L-μ3)
#     ∫_R = -(a1 * dnorm(R-μ1) + a2 * dnorm(R-μ2) + a3 * dnorm(R-μ3)) +
#           a1 * μ1 * pnorm(R-μ1) + a2 * μ2 * pnorm(R-μ2) + a3 * μ3 * pnorm(R-μ3)

#     μ_LR = (∫_R - ∫_L)/(pR - pL)

#     (pL, pR, ∫_L, ∫_R, μ_LR)
# end

"""
    cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)

Conditional mean on (L, R), defined as μ_{L,R} = (∫_L^R f(x) dx)^(-1) / (∫_L^R x f(x) dx)
"""
function cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)
    # (∫_l^r f(x) dx) = pR - pL
    pL = a1 * pnorm(L.-μ1) + a2 * pnorm(L.-μ2) #+ a3 * pnorm(L.-μ3)
    pR = a1 * pnorm(R.-μ1) + a2 * pnorm(R.-μ2) #+ a3 * pnorm(R.-μ3)

    # (∫_l^r x f(x) dx) = ∫_R - ∫_L
    ∫_L = -(a1 * dnorm(L.-μ1) + a2 * dnorm(L.-μ2)) + # + a3 * dnorm(L.-μ3)) +
          a1 * μ1 * pnorm(L.-μ1) + a2 * μ2 * pnorm(L.-μ2) #+ a3 * μ3 * pnorm(L.-μ3)
    ∫_R = -(a1 * dnorm(R-μ1) + a2 * dnorm(R-μ2)) + # + a3 * dnorm(R-μ3)) +
          a1 * μ1 * pnorm(R-μ1) + a2 * μ2 * pnorm(R-μ2) #+ a3 * μ3 * pnorm(R-μ3)
    ## adding a3/μ3 terms yields inaccurate result.

    # population conditional mean on (L, R)
    μ_LR = (∫_R .- ∫_L) ./ (pR .- pL)

    (pL, pR, ∫_L, ∫_R, μ_LR)
end

"""
    find_split(d, a1, a2, μ1, μ2, L, R;
               a3=0., μ3=0., find_deltas=false, find_split=!find_deltas)

Split the cluster (L, R) into subclusters that maximizes G_{L,R}.
Returns centered conditional mean (deltas) of the left and right subcluster
if `find_deltas` is true. Otherwise return the split point.
"""
function find_split(d, a1, a2, μ1, μ2, L, R;
                    a3=0., μ3=0., find_deltas=false, find_split=!find_deltas)
    # p, int, int_en evaluated at `d`
    ∫_d = -(a1 * dnorm(d.-μ1) + a2 * dnorm(d.-μ2)) + # + a3 * dnorm(d.-μ3)) +
          a1 * μ1 * pnorm(d.-μ1) + a2 * μ2 * pnorm(d.-μ2) #+ a3 * μ3 * pnorm(d.-μ3)
    p_d = a1 * pnorm(d.-μ1) + a2 * pnorm(d.-μ2) #+ a3 * pnorm(d.-μ3)

    pL, pR, ∫_L, ∫_R, μ_LR = cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)

    # find split point with normlized G'_{L,R} (=h)(supp. p.32)
    h = μ_LR .* (p_d .- pL).^2 + (pL + pR .- 2p_d) .* (∫_d .- ∫_L) - d .* (p_d .- pL) .* (pR .- p_d)
    split = d[2:end][h[2:end] .* h[1:end-1] .< 0]

    if find_deltas && (sum(μ1 .< split .< μ2) == 0)
        return(1., 1.)
    end

    # split point
    s = max(split[μ1 .< split .< μ2]...)

    if find_split
        return(s)
    else
        _, _, _, _, δ_1 = cond_mean_on_LR(L, s, a1, a2, a3, μ1, μ2, μ3)
        δ_1 -= (L + s)/2  # centering

        _, _, _, _, δ_2 = cond_mean_on_LR(s, R, a1, a2, a3, μ1, μ2, μ3)
        δ_2 -= (s + R)/2  # centering

        return(δ_1, δ_2)
    end

end

"""
    find_truncation(a1, a2, μ1, μ2,
                    R_ub=(μ2 + 8), L_lb=(μ1 - 8),
                    tuning_par1=0.01, tuning_par2=0.005,
                    a3=0, μ3=0, σ=1)

Find population right-truncation points for bimodal Normal mixtures.
Inputs:
    a1, a2, a3: weights for Normal distributions. a1 + a2 + a3 should be == 1
    μ1, μ2, μ3: means for normals μ1 = -μ2

Possible returns:
    0: negligible result (not a symmetric problem)
    -1: no truncation
"""
function find_truncation(a2, μ2;
                         a1=1-a2, μ1=-μ2, R_ub=μ2+8., L_lb=μ1-8.,
                         tuning_par1=0.01, tuning_par2=0.005,
                         a3=0., μ3=0., σ=1.)

    r = Array(R_ub:(-tuning_par1):μ1)

    # if the bimodal distn is symmetric, and a1 == a2,
    # then truncate both ends equally
    if a1 == a2 && μ2 == -μ1
        for i=1:length(r)

            # truncate both left and right with the same amount
            R = r[i]
            L = -R

            d = Array(L:0.01:R)
            δ_1, δ_2 = find_split(d, a1, a2, μ1, μ2, L, R; find_deltas=true)

            # L^* is the smallest L that satisfies δ_1 < 0 < δ_2. (supp. p.30)
            # Thus returns right-trunctation point
            if δ_1 < 0 < δ_2
                return(R)
            end
        end
    end

    # possible cases:
    #   a1 != a2 && μ2 != -μ1
    #   a1 != a2 && μ2 == -μ1
    if a1 != a2
        L_old = -1000
        for i=1:length(r)
            R = r[i]
            L = Array{Float64,1}(L_old:tuning_par2:(R-0.05))

            pL, pR, ∫_L, ∫_R, μ_LR = cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)

            # locate R_L using the equation μ_{L,R_L} - (L + R_L)/2.  (supp. p.30)
            temp = sign.(μ_LR .- (L .+ R)./2)

            if sum(temp[2:end] .* temp[1:end-1] .< 0) ≠ 0
                # this is L^* from (supp. p.28)
                L = L[2:end][temp[2:end] .* temp[1:end-1] .< 0][1]

                d = Array(L:0.01:R)
                δ_1, δ_2 = find_split(d, a1, a2, μ1, μ2, L, R; find_deltas=true)

                pL = a1 * pnorm(L-μ1) + a2 * pnorm(L-μ2) + a3 * pnorm(L-μ3)
                pR = a1 * pnorm(R-μ1) + a2 * pnorm(R-μ2) + a3 * pnorm(R-μ3)

                if sign((R-L)/(pR-pL) * (a1 * dnorm(L-μ1) + a2 * dnorm(L-μ2) + a3 * dnorm(L-μ3)) - 1) ≠ -1
                    println("🤔🤔 input is not good")
                end

                # L^* is the smallest L that satisfies δ_1 < 0 < δ_2. (supp. p.30)
                # Thus returns right-trunctation point
                if δ_1 < 0 < δ_2 && L_old ≤ L

                    # if L^* is achieved at i==1, it means there is no truncation.
                    if i == 1
                        return(-1)
                    end

                    return(R)
                end

                # increment L
                if L_old ≤ L
                    L_old = L
                end

                # if L > 0, it means roughly that L > m_ (m_≤0)  (supp. p.29 table and p.28)
                if L > 0
                    return(-1)
                end

            end
        end
    end

    # If a1 == a2 and μ1 != -μ2, then return 0
    # : do not need to assess in this case, since it can WLOG be reduced to a symmetric case)
    return(0)
end

"""
    clusterpath_pop(a2::T, μ2::T; a1=1-a2, μ1=-μ1, a3=0., μ3=0, σ=1) where T <: Real

Inputs:
    a2: weight on 2nd normal distribution
    μ2: mean of 2nd normal distribution
Returns:
    A dictionary of weights (`a1`, `a2`), means (`μ1`, `μ2`) and
    description of truncation (`L*`, `R*`) and split point (`s`).
"""
function clusterpath_pop(a2::T1, μ2::T2; a1=1-a2, μ1=-μ2, a3=0., μ3=0, σ=1.) where T1 <: Real where T2 <: Real
    # search right truncation point
    R = find_truncation(a2, μ2; a1=a1, μ1=μ1, a3=a3, μ3=μ3, σ=σ)

    # R == 0 or R == -1
    if R == 0 || R == -1
        d = Array(μ1:1e-04:μ2)
        # d[argmin(abs.(a1 * dnorm(d.-μ1) .- a2 * dnorm(d.-μ2)))]  ## ???????
        # d[argmin(a1 * dnorm(d.-μ1) .+ a2 * dnorm(d.-μ2))]

        return(Dict("a1" => a1, "a2" => a2,
                    "mu1" => μ1, "mu2" => μ2,
                    "s" => NaN, "L*" => NaN, "R*" => NaN))
    end

    # If R is a proper right-truncation point, then find split point
    L_lb = -3 * max(R, 3)  # lower bound of searching grid of L
    L_ub = R - 0.01  # upper bound
    L = Array(L_lb:5e-04:L_ub)

    _, _, _, _, μ_LR = cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)
    δs = sign.(μ_LR .- (L .+ R)/2)

    # L^* is the smallest L that satisfies δ_1 < 0 < δ_2
    L = L[2:end][δs[2:end] .* δs[1:end-1] .< 0][1]

    d = Array(L:0.001:R)
    split = find_split(d, a1, a2, μ1, μ2, L, R)

    return(Dict("a1" => a1, "a2" => a2,
                "mu1" => μ1, "mu2" => μ2,
                "s" => split, "L*" => L, "R*" => R))
end

# ---- simulation data generation ---- #

"""
    generate_mixture_normal(n, m, p)

generate `n` observations from mixture of univariate normals
each with standard deviation 1 and mean parameters `m` and proportion `p`.
"""
function generate_mixture_normal(n::Int64, m::Array{T1,1}=[-20., 0., 20.], p::Array{T2,1}=ones(length(m))/length(m)) where T1 <: Real where T2 <: Real
    @assert sum(p) ≈ 1
    mixed_normal = Array{Float64}(undef, 0)  # placeholder for the resulting mixture
    n_mixtures = floor.(Int64, n * p)  # sample sizes for each univariate normals

    # leftover n's is assigned to the first distribution > sum(n_mixtures)
    n_mixtures[1] += (n - sum(n_mixtures))

    # sample from Normal(0, 1)
    for i=1:length(m)
        mixed_normal = [mixed_normal; randn(n_mixtures[i]) .+ m[i]]
    end

    mixed_normal
end

end  # end module


text:
using Test

modname(fname)=splitext(basename(fname))[1]

#
# Include all Julia files in `testdir` whose name starts with `prefix`,
# Each file `prefixModName.jl` must contain a module named
# `prefixModName` which has a method test() returning true
# or false depending on success.
#
function run_tests_from_directory(testdir,prefix)
    println("Directory $(testdir):")
    @time begin
        examples=modname.(readdir(testdir))
        for example in examples
            println("  $(example):")
            if length(example)>=length(prefix) &&example[1:length(prefix)]==prefix
                path=joinpath(testdir,"$(example).jl")
                @eval begin
                    include($path)
                    # Compile + run test
                    print("   compile:")
                    @time @test eval(Meta.parse("$($example).test()"))
                    # Second run: pure execution time.
                    print("       run:")
                    @time eval(Meta.parse("$($example).test()"))
                end
            end
        end
    end
end

function run_all_tests()
    @time begin
        run_tests_from_directory(@__DIR__,"test_")
        run_tests_from_directory(joinpath(@__DIR__,"..","examples"),"Example")
    end
end

run_all_tests()


text:
export process_lifecycles

function collect_lifecycles(data::HONData, only_closed::Bool=false)
    simplices = data.simplices
    nverts = data.nverts
    times = data.times

    A, At, B = basic_matrices(simplices, nverts)
    simplex_order = simplex_degree_order(At)
    triangle_order = proj_graph_degree_order(B)

    all_events = Int64[]
    all_times = Int64[]
    num_events = Int64[]
    I, J, K = Int64[], Int64[], Int64[]

    for i = 1:size(B, 2)
        for (j, k) in combinations(neighbors(B, triangle_order, i), 2)
            if B[j, k] > 0
                # We have a triangle -- get the trail
                trail_events = Int64[]
                trail_times = Int64[]
                trail_simplices = Int64[]

                # SKIP if triangle is open and we only want to look at
                # closed triangles
                if only_closed && !triangle_closed(A, At, simplex_order, i, j, k)
                    continue
                end

                # Sort least common to most common
                a, b, c = i, j, k
                minval = min(simplex_order[i], simplex_order[j], simplex_order[k])
                if     minval == simplex_order[j]; a, b, c = j, i, k
                elseif minval == simplex_order[k]; a, b, c = k, i, j
                end

                # Get all simplices with a as an author
                for simplex_id in nz_row_inds(At, a)
                    # Check if b and c are in the simplex
                    has_b = A[b, simplex_id] > 0
                    has_c = A[c, simplex_id] > 0
                    if      has_b &&  has_c
                        push!(trail_events, 4)
                        push!(trail_times, times[simplex_id])
                        push!(trail_simplices, simplex_id)
                    elseif !has_b &&  has_c
                        push!(trail_events, 2)
                        push!(trail_times, times[simplex_id])
                        push!(trail_simplices, simplex_id)
                    elseif  has_b && !has_c
                        push!(trail_events, 1)
                        push!(trail_times, times[simplex_id])
                        push!(trail_simplices, simplex_id)
                    end
                end

                # Get all (b, c) simplices without a
                for simplex_id in nz_row_inds(At, b)
                    if A[c, simplex_id] > 0 && A[a, simplex_id] == 0
                        push!(trail_events, 3)
                        push!(trail_times, times[simplex_id])
                        push!(trail_simplices, simplex_id)
                    end
                end

                sp = sortperm(trail_times)
                trail_events = trail_events[sp]
                trail_times = trail_times[sp]
                trail_simplices = trail_simplices[sp]

                closure_ind = findnext(x -> x == 4, trail_events, 1)
                if closure_ind != nothing
                    trail_events = trail_events[1:closure_ind]
                    trail_times = trail_times[1:closure_ind]
                    trail_simplices = trail_simplices[1:closure_ind]
                end

                append!(all_events, trail_events)
                append!(all_times, trail_times)
                push!(num_events, length(trail_times))
                push!(I, a)
                push!(J, b)
                push!(K, c)
            end
        end
    end

    return (I, J, K, num_events, all_events, all_times)
end

function process_trail(events::Vector{Int64})
    ij, jk, ik = 0, 0, 0
    states = Int64[]
    push!(states, 1)

    for event in events
        curr_state = states[end]
        if     event == 1; ij += 1
        elseif event == 2; ik += 1
        elseif event == 3; jk += 1
        end

        # Process transition
        if event == 4
            push!(states, 11)  # filled triangle
            break
        end

        vals = [ij, ik, jk]
        sort!(vals)
        if     vals[2] == 0  # single edge
            if     vals[3] == 1; push!(states, 2)  # (0, 0, 1)
            else                 push!(states, 3)  # (0, 0, 2+)
            end
        elseif vals[1] == 0 && vals[2] > 0  # open wedge
            if     vals[2] == 1 && vals[3] == 1; push!(states, 4)  # (0, 1, 1)
            elseif vals[2] == 1 && vals[3] > 1;  push!(states, 5)  # (0, 1, 2+)
            elseif vals[2] > 1  && vals[3] > 1;  push!(states, 7)  # (0, 2+, 2+)
            else   error("Unknown transition")
            end
        elseif vals[1] > 0  # open triangle
            if     vals[1] == vals[2] == vals[3] == 1; push!(states, 6)  # (1, 1, 1)
            elseif vals[2] == 1 && vals[3] > 1;        push!(states, 8)  # (1, 1, 2+)
            elseif vals[1] == 1 && vals[2] > 1;        push!(states, 9)  # (1, 2+, 2+)
            elseif vals[1] > 1;                        push!(states, 10) # (2+, 2+, 2+)
            else error("Unknwon transition")
            end
        else
            error("Unknwon transition")
        end
    end
    return states
end

function earliest_activity(At::SpIntMat, times::Vector{Int64})
    earliest = zeros(Int64, size(At, 2))
    for j in 1:size(At, 2), time in times[nz_row_inds(At, j)]
        if earliest[j] == 0 || time < earliest[j]
            earliest[j] = time
        end
    end
    return earliest
end

"""
process_lifecycles
------------------

Process all lifecycles in a given dataset.

process_lifecycles(data::HONData)

Input parameters:
- data::HONData: The dataset

Returns a tuple (closed_transitions, open_transitions)
- closed_transitions::Array{Int64,2}: matrix whose (i, j) entry is the number of transitions from configuration j to configuration i, out of all triples of nodes that eventually simplicially close
- open_transitions::Array{Int64,2}: matrix whose (i, j) entry is the number of transitions from configuration j to configuration i, out of all triples of nodes that do not simplicially close
"""
function process_lifecycles(data::HONData)
    (triangles1, triangles2, triangles3, num_events, all_events, all_times) =
        collect_lifecycles(data, false)

    # Get earliest simplex of each node
    A, At, B = basic_matrices(data)
    earliest_times = earliest_activity(At, data.times)
    open_transitions = zeros(Int64, 11, 11)
    open_times = zeros(Int64, 11, 11)
    closed_transitions = zeros(Int64, 11, 11)
    closed_times = zeros(Int64, 11, 11)

    curr_ind = 1
    for (i, num) in enumerate(num_events)
        curr_events = all_events[curr_ind:(curr_ind + num - 1)]
        curr_times  = all_times[curr_ind:(curr_ind + num - 1)]
        curr_ind += num

        states = process_trail(curr_events)

        a = triangles1[i]
        b = triangles2[i]
        c = triangles3[i]
        # start time at the latest of the earliest simplex times
        emaxtime = maximum(earliest_times[[a, b, c]])
        if minimum(emaxtime) == 0
            error("0 Time (is value missing?)")
        end
        time = emaxtime

        for i = 2:length(states)
            s1 = states[i]
            s2 = states[i - 1]
            if states[end] == 11; closed_transitions[s1, s2] += 1
            else                  open_transitions[s1, s2]   += 1
            end
            if s2 != s1
                diff = curr_times[i - 1] - time
                if (s2 == 2 || s2 == 3) && s1 == 11  # single edge --> closed
                    diff = curr_times[i - 1] - max(emaxtime, time)
                end
                if (s1 == 2 && s2 == 1)  # empty --> single edge
                    diff = curr_times[i - 1] - emaxtime
                    if diff < 0
                        # Take the time since first co-authorship of 2 authors on the edge
                        next_start = median([earliest_times[a], earliest_times[b], earliest_times[c]])
                        diff = curr_times[i - 1] - next_start
                    end
                end
                if states[end] == 11; closed_times[s1, s2] += diff
                else                  open_times[s1, s2]   += diff
                end
                # Update time
                time = curr_times[i - 1]
            end
        end
    end
    return (closed_transitions, open_transitions)
end


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

JLLWrappers.@generate_main_file_header("FFMPEG")
JLLWrappers.@generate_main_file("FFMPEG", UUID("b22a6f82-2f65-5046-a5b2-351ab43fb4e5"))
end  # module FFMPEG_jll