text:
using MAT
using Random
include("../src/HyperLocal.jl")
include("../include/FlowSeed.jl")

## Read in the matrix, if it isn't read in already

s = time()
@time M = matread("../data/stackoverflow_answer_H.mat")
LabelMatrix = M["LabelMatrix"]
LabelNames = M["LabelNames"]
MainLabels = M["MainLabels"]
H = M["H"]
order = round.(Int64,vec(sum(H,dims=2)))
d = vec(sum(H,dims=1))
volA = sum(d)
m,n = size(H)
Ht = sparse(H')
toc = time()-s
println("Done loading things into memory in $toc seconds.")

##
Tconds = Vector{Float64}()
Tsizes = Vector{Int64}()
labels = Vector{Int64}()
for lab = 1:length(MainLabels)
    label = MainLabels[lab]
    T = findnz(LabelMatrix[:,label])[1]
    nT = length(T)
    condT, volT, cutT = tl_cond(H,T,d,1.0,volA,order)
    if condT < .2
        println("$nT \t $condT \t"*LabelNames[label])
        push!(Tconds,condT)
        push!(labels,label)
        push!(Tsizes,nT)
    end
end

##
lnum = length(labels)

# See outer parameters
ntimes = 1
epsilon = 1.0
delta = 5000.0
# seednum = 100
# grownum = 10000

# Output from HyperLocal, large delta
hl_pr = zeros(lnum,ntimes)
hl_re = zeros(lnum,ntimes)
hl_f1 = zeros(lnum,ntimes)
hl_time = zeros(lnum,ntimes)
hl_size = zeros(lnum,ntimes)
newS = zeros(lnum,ntimes)
hl_cond = zeros(lnum,ntimes)

# Output from HyperLocal, delta = 1.0
hl_pr1 = zeros(lnum,ntimes)
hl_re1 = zeros(lnum,ntimes)
hl_f11 = zeros(lnum,ntimes)
hl_time1 = zeros(lnum,ntimes)
hl_size1 = zeros(lnum,ntimes)
newS1 = zeros(lnum,ntimes)
hl_cond1 = zeros(lnum,ntimes)

# Output from first baseline 1
b1_pr = zeros(lnum,ntimes)
b1_re = zeros(lnum,ntimes)
b1_f1 = zeros(lnum,ntimes)
b1_cond = zeros(lnum,ntimes)

# Output from baseline 2
b2_pr = zeros(lnum,ntimes)
b2_re = zeros(lnum,ntimes)
b2_f1 = zeros(lnum,ntimes)
b2_cond = zeros(lnum,ntimes)

# Keep track of R
r_pr =  zeros(lnum,ntimes)
r_re =  zeros(lnum,ntimes)
r_f1 =  zeros(lnum,ntimes)
r_cond =  zeros(lnum,ntimes)

sl_pr = zeros(lnum,ntimes)
sl_re = zeros(lnum,ntimes)
sl_f1 =zeros(lnum,ntimes)
sl_cond = zeros(lnum,ntimes)
sl_size = zeros(lnum,ntimes)
sl_time = zeros(lnum,ntimes)

zl_pr = zeros(lnum,ntimes)
zl_re = zeros(lnum,ntimes)
zl_f1 = zeros(lnum,ntimes)
zl_cond = zeros(lnum,ntimes)
zl_size = zeros(lnum,ntimes)
zl_time = zeros(lnum,ntimes)

# For each epsilon we store a different matrix of outputs
outputmat = "Output_Stack/Set45_$(epsilon)_$(delta).mat"
println(outputmat)

maz = matread("ZCE_Stack.mat")
Az = maz["Az"]
mas = matread("SCE_Stack.mat")
As = mas["As"]

## For a fixed epsilon, seednum, and grownum, run experiments
# on each cluster multiple times

for lab = 1:length(labels)
    label = labels[lab]
    T = findnz(LabelMatrix[:,label])[1]
    nT = length(T)
    @show nT,label

    for index = 1:ntimes

        # Generate a new seed set
        seednum = round(Int64,nT/20)
        grownum = round(Int64,min(nT*2,n))
        p = randperm(nT)
        Rstart = T[p[1:seednum]]
        OneHop = get_immediate_neighbors(H,Ht,Rstart)
        Rmore = BestNeighbors(H,d,Rstart,OneHop,grownum)
        R = union(Rmore,Rstart)
        Rs = findall(x->in(x,Rstart),R)     # Force seed nodes to be in output set
        prr, rer, f1r = PRF(T,R)
        r_pr[lab,index] = prr
        r_re[lab,index] = rer
        r_f1[lab,index] = f1r
        condR, volR, cutR = tl_cond(H,R,d,1.0,volA,order)
        r_cond[lab,index] = condR

        # Run HyperLocal with delta = 1.0
        s = time()
        S, lcond = HyperLocal(H,Ht,order,d,R,epsilon,1.0,Rs,true)
        hl_time1[lab,index] = time()-s
        condS, volS, cutS = tl_cond(H,S,d,1.0,volA,order)
        pr, re, f1_d1 = PRF(T,S)
        hl_pr1[lab,index] = pr
        hl_re1[lab,index] = re
        hl_f11[lab,index] = f1_d1
        hl_size1[lab,index] = length(S)
        hl_cond1[lab,index] = condS
        nS =  length(setdiff(S,R))
        newS1[lab,index] = nS

        # Run HyperLocal with delta = 1000
        s = time()
        S, lcond = HyperLocal(H,Ht,order,d,R,epsilon,delta,Rs,true)
        hl_time[lab,index] = time()-s
        condS, volS, cutS = tl_cond(H,S,d,1.0,volA,order)
        pr, re, f1 = PRF(T,S)
        hl_pr[lab,index] = pr
        hl_re[lab,index] = re
        hl_f1[lab,index] = f1
        hl_size[lab,index] = length(S)
        hl_cond[lab,index] = condS
        nS =  length(setdiff(S,R))
        newS[lab,index] = nS

        # First baseline
        kS = nT-length(Rstart)
        B1 = BestNeighbors(H,d,Rstart,OneHop,kS)
        pr1, re1, f11 = PRF(T,B1)
        b1_pr[lab,index] = pr1
        b1_re[lab,index] = re1
        b1_f1[lab,index] = f11

        # Baseline 2
        B2 = TopNeighbors(H,Rstart,OneHop,kS)
        pr2, re2, f12 = PRF(T,B2)
        b2_pr[lab,index] = pr2
        b2_re[lab,index] = re2
        b2_f1[lab,index] = f12

        # Simple Clique Expansion
        nR = length(R)
        Rs_vec = zeros(length(R))
        Rs_vec[Rs] .= 1
        starter = time()
        SL, lcond = FlowSeed(As,R,epsilon,zeros(nR),Rs_vec)
        sl_time[lab,index] = time()-starter
        pr3, re3, f13 = PRF(T,SL)
        condS, volS, cutS = tl_cond(H,SL,d,1.0,volA,order)
        sl_pr[lab,index] = pr3
        sl_re[lab,index] = re3
        sl_f1[lab,index] = f13
        sl_size[lab,index] = length(SL)
        sl_cond[lab,index] = condS

        # Zhou Clique Expansion + SimpleLocal
        starter = time()
        ZL, lcond = FlowSeed(Az,R,epsilon,zeros(nR),Rs_vec)
        zl_time[lab,index] = time()-starter
        pr4, re4, f14 = PRF(T,ZL)
        condS, volS, cutS = tl_cond(H,ZL,d,1.0,volA,order)
        zl_pr[lab,index] = pr4
        zl_re[lab,index] = re4
        zl_f1[lab,index] = f14
        zl_size[lab,index] = length(ZL)
        zl_cond[lab,index] = condS

        println("$label ($nT):\n $f1r \n $f11 \n $f12 \n $f1 \n $f1_d1\n \t $nS "*LabelNames[label])
    end

    matwrite(outputmat, Dict("hl_size"=>hl_size, "newS"=>newS, "hl_time"=>hl_time,
    "hl_pr"=>hl_pr, "hl_re"=>hl_re, "hl_f1"=>hl_f1, "hl_cond"=>hl_cond,
    "hl_pr1"=>hl_pr1, "hl_re1"=>hl_re1, "hl_f11"=>hl_f11, "hl_cond1"=>hl_cond1,
    "sl_pr"=>sl_pr, "sl_re"=>sl_re, "sl_f1"=>sl_f1, "sl_cond"=>sl_cond, "sl_time"=>sl_time,
    "zl_pr"=>zl_pr, "zl_re"=>zl_re, "zl_f1"=>zl_f1, "zl_cond"=>zl_cond, "zl_time"=>zl_time,
    "sl_size"=>sl_size, "zl_size"=>zl_size,
    "b1_pr"=>b1_pr, "b1_re"=>b1_re, "b1_f1"=>b1_f1,
    "r_pr"=>r_pr, "r_re"=>r_re, "r_f1"=>r_f1,
    "b2_pr"=>b2_pr, "b2_re"=>b2_re, "b2_f1"=>b2_f1))

end

max_stars_repo_name: nveldt/HypergraphFlowClustering

max_issues_repo_name: nveldt/HypergraphFlowClustering

max_forks_repo_name: nveldt/HypergraphFlowClustering


text:
##########################
# FunctionCell Interface #
##########################
"""
	libsvm([d] [;kwargs])

Generates a function cell that when piped data, trains a SVM classifier/regressor using the LIBSVM.jl interface.
`d::Distances.PreMetric` is a distance function that can be used to pre-compute the kernel.

# Keyword arguments (from LIBSVM.jl, defaults modified in some cases)
  * `svmtype::Type=LIBSVM.SVC`: Type of SVM to train `SVC` (for C-SVM), `NuSVC`,
    `OneClassSVM`, `EpsilonSVR` or `NuSVR`. Defaults to `OneClassSVM` if
      no labels are provided.
  * `kernel::Kernel.KERNEL=Kernel.Linear`; available: `Linear`, `Polynomial`
    `RadialBasis`, `Sigmoid` or `Precomputed`. If `d` is provided, the default is `Kernel.Precomputed`
  * `degree::Integer=1`: Kernel degree. Used for polynomial kernel
  * `gamma::Float64=1.0/size(X, 1)` : γ for kernels
  * `coef0::Float64=0.0`: parameter for sigmoid and polynomial kernel
  * `cost::Float64=1.0`: cost parameter C of C-SVC, epsilon-SVR, and nu-SVR
  * `nu::Float64=0.5`: parameter nu of nu-SVC, one-class SVM, and nu-SVR
  * `epsilon::Float64=0.1`: epsilon in loss function of epsilon-SVR
  * `tolerance::Float64=0.001`: tolerance of termination criterion
  * `shrinking::Bool=true`: whether to use the shrinking heuristics
  * `probability::Bool=true`: whether to train a SVC or SVR model for probability estimates
  * `weights::Union{Dict{T, Float64}, Void}=nothing`: dictionary of class weights
  * `cachesize::Float64=4000.0`: cache memory size in MB
  * `verbose::Bool=false`: print training output from LIBSVM if true

!!! note
	For regression (`svmtype=LIBSVM.NuSVR`, or `svmtype=LIBSVM.EpsilonSVR`) `probability`
	must be set explicitly to `false` otherwise Julia might crash (e.g. unsafe_copy! problem).
	The problem is present in Julia 0.6 and LIBSVM v.0.1.0

Read the `LIBSVM.jl` documentation for more information.
"""
libsvm(;kwargs...) = FunctionCell(libsvm, (), ModelProperties(), kwtitle("LIBSVM",kwargs); _libsvm_parse_kwargs_(nothing, kwargs)...)
libsvm(d::T where T<:Distances.PreMetric; kwargs...) = FunctionCell(libsvm, (d,), ModelProperties(),  kwtitle("LIBSVM",kwargs); _libsvm_parse_kwargs_(d, kwargs)...)

############################
# DataCell/Array Interface #
############################
"""
	libsvm(x, [,d] [;kwargs])

Generates a trained function cell by training a SVM classifier/regressor using the LIBSVM.jl interface.
`x`::DataCell is the training data and `d::Distances.PreMetric` is a distance function that can be used
to pre-compute the kernel.
"""
# Training
libsvm(x::T where T<:CellData; kwargs...) = libsvm((getx!(x), gety(x)), nothing; _libsvm_parse_kwargs_(nothing, kwargs)...)

libsvm(x::T where T<:CellData, d::S where S<:Distances.PreMetric; kwargs...) = libsvm((getx!(x), gety(x)), d; _libsvm_parse_kwargs_(d, kwargs)...)

libsvm(x::T where T<:AbstractArray, d::Union{Void, Distances.PreMetric}=nothing; kwargs...) = libsvm((x, nothing), d; _libsvm_parse_kwargs_(d, kwargs)...)

libsvm(x::Tuple{T,S} where T<:AbstractVector where S<:Union{Void,AbstractVector}, d::Union{Void, Distances.PreMetric}=nothing; kwargs...) =
	libsvm((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), d; _libsvm_parse_kwargs_(d, kwargs)...)

libsvm(x::Tuple{T,S} where T<:AbstractMatrix where S<:Union{Void,AbstractVector}, d::Union{Void, Distances.PreMetric}=nothing; kwargs...) = begin

	newkwargs = _libsvm_parse_kwargs_(d, kwargs)

	# If no labels present, change to OneClassSVM, otherwise check label size
	if x[2] isa Void
		newkwargs[:svmtype] = LIBSVM.OneClassSVM
	else
		@assert nobs(x[1]) == nobs(x[2]) "[libsvm] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."
	end

	# Transform labels first if svmtype if SVC or NuSVC
	yenc=[];
	clf = true
	if isequal(newkwargs[:svmtype], LIBSVM.SVC) || isequal(newkwargs[:svmtype], LIBSVM.NuSVC)
		enc = labelencn(x[2])
		yenc = label2ind.(x[2],enc)
	else
		enc = labelencn([])
		yenc = x[2] # regression or oneclass svm, do not change
		clf = false
	end

	# Train model
	if d isa Void
		# Use internal kernel
		if isequal(newkwargs[:svmtype], LIBSVM.OneClassSVM)
			modeldata = LIBSVM.svmtrain(getobs(x[1]); newkwargs...) 					# one-class SVM
		else
			modeldata = LIBSVM.svmtrain(getobs(x[1]), yenc; newkwargs...)					# classifier/regressor
		end
	else
		# Use precomputed kernels
		if isequal(newkwargs[:svmtype], LIBSVM.OneClassSVM)
			modeldata = LIBSVM.svmtrain([collect(1:nobs(x[1]))';x[1] |> dist(x[1])]; newkwargs...) 		# one-class SVM
		else
			modeldata = LIBSVM.svmtrain([collect(1:nobs(x[1]))';x[1] |> dist(x[1])], yenc; newkwargs...) 	# classifier/regressor
		end
	end

	# Keep track of how the labeles were re-ordered
	reorder = indexin(1:length(enc.label), modeldata.labels) # it is possible to use 1:length(enc.label) because the unique labels are sorted
							  	 # and the labels are encoded according to the sorted unique labels (i.e. position)
	# Build model properties
	modelprops = ModelProperties(nvars(x[1]),
		   	  	     clf ? length(enc.label) : 1, # the output size is the number of classes (e.g. 1 prob/class) for classification, 1 for regression
			  	     enc,
				     Dict("kernel_distance" => d)
	)

	# Build title
	title = "LIBSVM: svmtype=$(newkwargs[:svmtype]), kernel=$(newkwargs[:kernel]), probability=$(newkwargs[:probability]), cachesize=$(newkwargs[:cachesize])(MB)"

	# Return trained cell (retain both SVM model and initial training data)
	FunctionCell(libsvm, Model((modeldata, getobs(x[1]), reorder), modelprops), title)
end

# Execution
libsvm(x::T where T<:CellData, model::Model{<:Tuple}) =
	datacell(libsvm(getx!(x), model), gety(x))
libsvm(x::T where T<:AbstractVector, model::Model{<:Tuple}) =
	libsvm(mat(x, LearnBase.ObsDim.Constant{2}()), model)
libsvm(x::T where T<:AbstractMatrix, model::Model{<:Tuple})::Matrix{Float64} =
	_libsvm_execute_(model.data[1].SVMtype, model.data[1], x, model.properties.other["kernel_distance"], model.data[2], model.data[3])

# Define execution methods (for return type annotations) and coversion of results to matrices
_libsvm_execute_(::Type{LIBSVM.OneClassSVM}, model_data, x, ::Void, xd, reorder=Int[]) =
	mat(LIBSVM.svmpredict(model_data, x)[1], LearnBase.ObsDim.Constant{2}())

_libsvm_execute_(::Type{LIBSVM.OneClassSVM}, model_data, x, d, xd, reorder=Int[]) =
	mat(LIBSVM.svmpredict(model_data, [collect(1:nobs(x))'; x |> dist(xd, d)])[1], LearnBase.ObsDim.Constant{2}())

_libsvm_execute_(::Union{Type{LIBSVM.EpsilonSVR}, Type{LIBSVM.NuSVR}}, model_data, x, ::Void, xdm, reorder=Int[]) =
	mat(LIBSVM.svmpredict(model_data, x)[1], LearnBase.ObsDim.Constant{2}())

_libsvm_execute_(::Union{Type{LIBSVM.EpsilonSVR}, Type{LIBSVM.NuSVR}}, model_data, x, d, xd, reorder=Int[]) =
	mat(LIBSVM.svmpredict(model_data, [collect(1:nobs(x))'; x |> dist(xd, d)])[1], LearnBase.ObsDim.Constant{2}())

_libsvm_execute_(::Union{Type{LIBSVM.SVC}, Type{LIBSVM.NuSVC}}, model_data, x, ::Void, xd, reorder::Vector{Int}) =
	LIBSVM.svmpredict(model_data, x)[2][reorder,:]

_libsvm_execute_(::Union{Type{LIBSVM.SVC}, Type{LIBSVM.NuSVC}}, model_data, x, d, xd, reorder::Vector{Int}) =
	LIBSVM.svmpredict(model_data, [collect(1:nobs(x))'; x |> dist(xd, d)])[2][reorder,:]

# A function that generates default keyword arguments
function _libsvm_parse_kwargs_(d, kwargs)

	outkwargs = Dict(
	   (:svmtype => LIBSVM.SVC), # SVM type
	   (:kernel => d isa Void ? LIBSVM.Kernel.Linear : LIBSVM.Kernel.Precomputed),
	   (:degree => 1),
	   #(:gamma => 1.0/size(X,1)),
	   (:coef0 => 0.0),
	   (:cost => 1.0),
	   (:nu => 0.5),
	   (:epsilon => 0.1),
	   (:tolerance => 0.001),
	   (:shrinking => true),
	   (:probability => true),
	   (:weights => nothing),
	   (:cachesize => 4000.0), 	# cache 4GBs
	   (:verbose => false)
	)

	# Loop through any keyword arguments and modify the defaults
	for (kwarg, argval) in kwargs
		if kwarg == :svmtype
			# Check if we have regression, if so , make probability=false
			if argval != LIBSVM.SVC && argval != LIBSVM.NuSVC
				outkwargs[:probability] = false
			end
		end

		# Update all default arguments with new values
		if kwarg in keys(outkwargs)
			outkwargs[kwarg] = argval
		else
			push!(outkwargs, kwarg=>argval)
		end
	end

	# Return concatenation of defaults and keyword arguments; any arguments from the
	# defaults that appear in kwargs explicitly will be modified, as the should be.
	return outkwargs
end

ext: jl


text:
cd(dirname(Base.source_path()))
push!(LOAD_PATH, dirname(Base.source_path()))

workspace()
reload("JuDaisy")
using JuDaisy

# im must be transposed, since julia stores it in column-major order
im = Float32[0.4902 0.4824 0.5137 0.5529 0.5451 0.5490 0.5490 0.5804 0.6039 0.6039 0.5882 0.5961 0.5176
             0.5020 0.4941 0.5255 0.5569 0.5529 0.5647 0.5647 0.5725 0.5882 0.6039 0.5882 0.6078 0.5804
             0.4980 0.5333 0.5882 0.5843 0.5882 0.5882 0.5765 0.5922 0.5882 0.5961 0.6118 0.6039 0.6000
             0.4980 0.5843 0.6078 0.5843 0.5608 0.6000 0.6000 0.5804 0.6078 0.5922 0.6118 0.6078 0.6000
             0.5647 0.5804 0.6118 0.5608 0.5647 0.5725 0.6157 0.5608 0.6000 0.6078 0.6157 0.6000 0.5922
             0.5686 0.5882 0.5961 0.6078 0.7216 0.7490 0.7216 0.7373 0.5961 0.6078 0.6196 0.6078 0.6039
             0.5569 0.5451 0.6078 0.6941 0.7098 0.8275 0.7608 0.7020 0.5804 0.5412 0.5882 0.6000 0.5529
             0.5922 0.5882 0.7922 0.6706 0.6941 0.8118 0.6275 0.6824 0.5176 0.4745 0.5725 0.5882 0.5765
             0.5412 0.5451 0.7255 0.6118 0.7765 0.7961 0.7137 0.7294 0.4510 0.4588 0.5882 0.5725 0.5843
             0.9882 0.8235 0.6235 0.6863 0.7490 0.6980 0.6314 0.4706 0.4824 0.4471 0.5882 0.5608 0.6118
             0.8588 0.9765 0.8627 0.6627 0.6824 0.6392 0.6275 0.5059 0.4314 0.5137 0.5451 0.5490 0.5490
             0.7373 0.9529 0.8431 0.5843 0.6000 0.5529 0.4667 0.4314 0.5098 0.5529 0.5373 0.5882 0.5843
             0.6863 0.8941 0.8941 0.5373 0.5490 0.5451 0.4510 0.4078 0.5098 0.4941 0.5529 0.5294 0.5412
             0.6902 0.7255 0.6588 0.4549 0.5020 0.5059 0.4980 0.3804 0.4980 0.5608 0.6275 0.5765 0.6000]';

using Images, ImageView
prop = {"colorspace" => "Gray", "spatialorder" => ["x", "y"]}
imA = Image(im, prop);
show(imA)
imgc, imgslice = view(imA, pixelspacing = [1,1]); Tk.set_size(toplevel(imgc), 300, 300)

desc = Daisy()

image!(desc, im)
parameters!(desc, 3.)
@show desc.descriptor_size

compute!(desc)
normalize!(desc)

display(desc[8,8])


text:
# Include packages
if !(@isdefined peaks)
    include("peaks.jl")
end
using Printf: @printf, @sprintf

####################
## function
####################
function grad2d(dataorg::Array{Float64,2})
    ny, nx = size(dataorg)
    Gx = zeros(ny,nx)
    Gy = zeros(ny,nx)

    diffx = diff(dataorg, dims=2)
    Gx[:,1], Gx[:,end] = diffx[:,1], diffx[:,end]
    Gx[:,2:nx-1] = [(dataorg[i,j+1]-dataorg[i,j-1])/2 for i=1:ny, j=2:nx-1]

    diffy = diff(dataorg, dims=1)
    Gy[1,:], Gy[end,:] = diffy[1,:], diffy[end,:]
    Gy[2:ny-1,:] = [(dataorg[i+1,j]-dataorg[i-1,j])/2 for i=2:ny-1, j=1:nx]

    return(Gx, Gy)
end
####################

####################
## main
####################

# Parameters
const N=31
qs = Int64(2)
L=0.25
# create mesh data
xmat, ymat, ϕ = peaks(N)
xvec = vec(xmat[1,:])
yvec = vec(ymat[:,1])
# gradient
ϕx,ϕy = grad2d(ϕ)
ϕxx,_ = grad2d(ϕx)
_,ϕyy = grad2d(ϕy)

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

####################
## plot
####################

# directory where figures are printed
figdir="./fig"
if !isdir(figdir); mkdir(figdir); end

using Plots
pyplot()
#gr()

########
## Figure 1
# Two-dimensional non-filled contour
clibrary(:misc)
contour(xvec, yvec, ϕ, c=(:rainbow), fill=false, tickfont=12,
        xlims=(-3.,3.), ylims=(-3.,3.), clims=(-6.,8.),
        xlabel="X", ylabel="Y", axis_ratio=:equal, size=(800,600),
        )
        # arrows
quiver!(vec(xmat[1:qs:end,1:qs:end]), vec(ymat[1:qs:end,1:qs:end]),
        quiver=(-L*vec(ϕx[1:qs:end,1:qs:end]), -L*vec(ϕy[1:qs:end,1:qs:end])),
        color=:black, arrow=arrow(:closed,:head,0.3,0.3))
# save
savefig(joinpath(figdir,"contour_grad2d.png"))
########

########
## Figure 2
# filled contour
clibrary(:misc)
contour(xvec, yvec, ϕ, c=(:rainbow), fill=true, tickfont=12,
        xlims=(-3.,3.), ylims=(-3.,3.), clims=(-6.,8.),
        xlabel="X", ylabel="Y", axis_ratio=:equal, size=(800,600),
        )
# arrows
quiver!(vec(xmat[1:qs:end,1:qs:end]), vec(ymat[1:qs:end,1:qs:end]),
        quiver=(L*vec(ϕxx[1:qs:end,1:qs:end]), L*vec(ϕyy[1:qs:end,1:qs:end])),
        color=:black, arrow=arrow(:closed,:head,0.3,0.3))
# save
savefig(joinpath(figdir,"contour-filled_grad2d.png"))
########

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


text:
@testset "Caching" begin
  # @testset "Empty string" begin
  #   using Genie, Genie.Cache

  #   @test String(r.body) == ""
  # end;

end;


text:
module MLMetrics_Test
using MLLabelUtils
using Test, ReferenceTests

# check for ambiguities
refambs = detect_ambiguities(Base, Core)
using MLMetrics
ambs = detect_ambiguities(MLMetrics, Base, Core)
@test Set(setdiff(ambs, refambs)) == Set{Tuple{Method,Method}}()

classification_tests = joinpath.("classification", [
    "counts.jl",
    "metrics.jl",
    "roc.jl",
])

for test in classification_tests
    @testset "$test" begin
        include(test)
    end
end

regression_tests = joinpath.("regression", [
    "regression.jl"
])

for test in regression_tests
    @testset "$test" begin
        include(test)
    end
end
end


text:
@testset "Inf" begin
  @test isinf(ArbFloat(Inf))
  @test isinf(ArbReal(Inf))
  @test isinf(ArbComplex(Inf))

  @test isposinf(ArbFloat(Inf))
  @test isposinf(ArbReal(Inf))
  @test isposinf(ArbComplex(Inf))

  @test isneginf(ArbFloat(-Inf))
  @test isneginf(ArbReal(-Inf))
  @test isneginf(ArbComplex(-Inf))
end

@testset "NaN" begin
  @test isnan(ArbFloat(NaN))
  @test isnan(ArbReal(NaN))
  @test isnan(ArbComplex(NaN))
end

@testset "01" begin
  @test iszero(ArbFloat(0))
  @test iszero(ArbReal(0))
  @test iszero(ArbComplex(0))

  @test isone(ArbFloat(1))
  @test isone(ArbReal(1))
  @test isone(ArbComplex(1))
end


text:
#****************************************************************#
# VeysiADN 4 Apr 2021
# Advanced Robotics Homework 3 Codes
#****************************************************************#
using LinearAlgebra
#****************************************************************#
# Returns 3 x 3 rotation matrix representing
# rotation about unit vector w by angle theta.
# theta should be in radian
function Rodrigues(w,theta)
    if( length(w) != 3)
        println("!Warning! : W matrix should have 3 elements.")
        return
    end
    if(size(w)==(1,3))
        w=w';
    end
    w_hat = [0      -w[3]    w[2];
             w[3]     0     -w[1];
            -w[2]    w[1]      0];
    rotation_matrix = cos(theta)*Matrix(1I,3,3) + sin(theta)*w_hat+(1-cos(theta))*w*w';
    return rotation_matrix;
end
#****************************************************************#
# Builds rigid transformation matrix T from given
#rotation matrix R and translation t
function RpToTrans(R,t)
    if( size(R)!=(3,3) )
        println("!Warning! : Rotation matrix should be 3x3.")
        return
    end
    if( length(t) != 3 )
        println("!Warning! : Translation matrix should have 3 elements.")
        return
    end
    if(size(t)==(1,3))
        t=t';
    end
    R_B_T = hcat(R,t) ;
    R_B_T = vcat(R_B_T,[0 0 0 1]) ;
    return R_B_T
end
#****************************************************************#
#Extracts rotation matrix R and translation t from given /
# rigid transformation matrix T
function TransToRp(T)
    if( size(T)!=(4,4) )
        println("!Warning! : Transformation matrix should be 4x4.")
        return
    end
    R = T[1:3,1:3]
    t = T[1:3,4]
    return R,t;
end
#****************************************************************#
# Inverts rigid transformation matrix T
function TransInv(T)
    if( size(T)!=(4,4) )
        println("!Warning! : W matrix should be 4x4.")
        return
    end
    if( det(T)==0 )
        println("!Warning! : Inverse doesn't exist (Determinant is equal to zero).")
        return
    end
    R,t = TransToRp(T);
    r_inv =hcat(R',-(R'*t))
    r_inv = vcat(r_inv,[0 0 0 1]);
    return r_inv;
end
#****************************************************************#

max_stars_repo_name: veysiadn/AdvancedRoboticsJulia

max_issues_repo_name: veysiadn/AdvancedRoboticsJulia

max_forks_repo_name: veysiadn/AdvancedRoboticsJulia


text:
using TropicalTensors
using Test

@testset "dumpload" begin
    include("dumpload.jl")
end

@testset "yao" begin
    include("yao.jl")
end

@testset "potts" begin
    include("potts.jl")
end

@testset "reversible" begin
    include("reversible/reversible.jl")
end

@testset "cuda" begin
    if Base.find_package("CuYao") !== nothing
        include("cuda.jl")
    end
end


text:
#=
## This is not needed anymore, because of the constructor redefinition to accept x::Integers and converts T(x) instead of accepting x::T.
This automatically enables the zero and one functions.

# creates some special elements
import Base.zero

"""
    zero(::ℤd{T, d})
or
    zero(x::ℤd{T, d})

Returns the zero element ℤd{T, d}(T(0)).
"""
zero(::ℤd{T, d}) where {T <: Integer, d} = ℤd{T, d}(T(0))
zero(x::ℤd{T, d}) where {T <: Integer, d} = ℤd{T, d}(T(0))
=#


text:
@testset "Lennard-Jones" begin
    # The parameters in this test come from the NIST computed values
    # in https://mmlapps.nist.gov/srs/LJ_PURE/md.htm
    rho = 0.84
    nrho = 200
    p = OrnsteinZernike.Parameters(1.0, 2.5, rho, 2^11, nrho, 0.85)
    pot = OrnsteinZernike.LennardJones()
    st = OrnsteinZernike.Structure(p, pot)
    cls = OrnsteinZernike.HypernettedChain()

    inter = OrnsteinZernike.Interaction(p, st, cls)

    result = OrnsteinZernike.solve(inter)
    if 1.0 ∈ result.r
        first_contact = result.gr[result.r .== 1.0][1]
    else
        first_contact = maximum(result.gr)
    end
    @show first_contact
    # We don't know exactly the value at contact, so just check that it is a valid number
    @test first_contact isa Real
end

function sw_oz(pot)
    rho = 0.8
    nrho = 30
    ktemp = 1.5
    p = OrnsteinZernike.Parameters(1.0, 6.0, rho, 2^12, nrho, ktemp)
    st = OrnsteinZernike.Structure(p, pot)
    cls = OrnsteinZernike.MeanSpherical()

    inter = OrnsteinZernike.Interaction(p, st, cls)

    result = OrnsteinZernike.solve(inter)

    return result
end

@testset "Square Well" begin
    potentials = Dict("1.5" => SquareWell(1.5), "1.3" => SquareWell(1.3))

    for (k, v) in potentials
        @testset "$k" begin
            result = sw_oz(v)
            if 1.0 ∈ result.r
                first_contact = result.gr[result.r .== 1.0][1]
            else
                first_contact = maximum(result.gr)
            end
            @show "SW $k" first_contact
            @test first_contact isa Real
        end
    end
end

function phs_oz(pot)
    rho = 0.8
    nrho = 30
    p = OrnsteinZernike.Parameters(1.0, 6.0, rho, 2^11, nrho, 1.0)
    st = OrnsteinZernike.Structure(p, pot)
    cls = OrnsteinZernike.ModifiedVerlet()

    inter = OrnsteinZernike.Interaction(p, st, cls)

    return OrnsteinZernike.solve(inter)
end

@testset "Pseudo Hard Sphere" begin
    potentials = Dict("50" => PseudoHS(), "51" => PseudoHS(51.0))

    for (k, v) in potentials
        @testset "$k" begin
            result = phs_oz(v)
            if 1.0 ∈ result.r
                first_contact = result.gr[result.r .== 1.0][1]
            else
                first_contact = maximum(result.gr)
            end
            @show "PHS $k" first_contact
            @test first_contact isa Real
        end
    end
end


text:
# Polyhedral gravitation

using LinearAlgebra
using StaticArrays
using HalfEdges

function setmesh(shapepath)
    topology, points = loadmesh(shapepath)
    # points *= 1000.

    println("Vertices : ", nvertices(topology))
    println("Faces    : ", nfaces(topology))
    println("Edges    : ", nedges(topology))
end

"""
    getFaceCenters(topology, points) -> facecenters

Get an array of center-positions of all faces

# Arguments
- `topology` : Topology instace of HalfEdges.jl
- `points`   : vertex-list of the polyhedron
# Return
- facecenters : array of center-positions of faces
"""
function getFaceCenters(topology, points)
    num_face = nfaces(topology)
    faces = facelist(topology)
    facecenters = Array{SVector{3,Float64}, 1}(undef, num_face)

    for i in 1:num_face
        face = faces[i]
        center = points[face[1]] + points[face[2]] + points[face[3]]
        center /= 3
        facecenters[i] = center
    end
    return facecenters
end

"""
    getFaceAreas(topology, points) -> faceareas

Get an array of areas of all faces

# Arguments
- `topology` : Topology instace of HalfEdges.jl
- `points`   : vertex-list of the polyhedron
# Return
- faceareas : array of areas of faces
"""
function getFaceAreas(topology, points)
    num_face = nfaces(topology)
    faceareas = Array{Float64, 1}(undef, num_face)

    for i in 1:num_face
        area = HalfEdges.area(topology, points, FaceHandle(i))
        faceareas[i] = area
    end
    return faceareas
end

"""
    getFaceNormals(topology, points) -> facenormals

Get an array of normal vectors to all faces

# Arguments
- `topology` : Topology instace of HalfEdges.jl
- `points`   : vertex-list of the polyhedron
# Return
- facenormals : array of face normals
"""
getFaceNormals(topology, points) = HalfEdges.normals(topology, points)

"""
    getFaceVertices(topology, points) -> r1s, r2s, r3s
"""
function getFaceVertices(topology, points)
    num_face = nfaces(topology)
    faces = facelist(topology)

    r1s = Array{SVector{3,Float64}, 1}(undef, num_face)
    r2s = Array{SVector{3,Float64}, 1}(undef, num_face)
    r3s = Array{SVector{3,Float64}, 1}(undef, num_face)

    for i in 1:num_face
        r1s[i] = points[faces[i][1]]
        r2s[i] = points[faces[i][2]]
        r3s[i] = points[faces[i][3]]
    end
    return r1s, r2s, r3s
end

"""
    getHalfEdgeVertices(topology, points) -> heads, tails
"""
function getHalfEdgeVertices(topology, points)
    num_halfedge = nhalfedges(topology)

    heads = Array{SVector{3,Float64}, 1}(undef, num_halfedge)
    tails = Array{SVector{3,Float64}, 1}(undef, num_halfedge)

    for i in 1:num_halfedge
        he = topology.he[i]

        heads[i] = points[he.head]                       # head of the halfedge
        tails[i] = points[HalfEdges.tail(topology, he)]  # tail of the halfedge
    end
    return heads, tails
end

################################################################
#                Constant density polyhedron
################################################################

getSurfaceGravity(topology, points, ρ) = getgravity(topology, points, ρ, getFaceCenters(topology, points))

function getgravity(topology, points, ρ, positions)
    num_halfedge = nhalfedges(topology)
    num_face = nfaces(topology)
    num_data = length(positions)

    potentials = zeros(num_data)
    accelerations = [SA_F64[0, 0, 0] for i in 1:num_data]
    # accelerations = [[0., 0., 0.] for i in 1:num_data]

    heads, tails = getHalfEdgeVertices(topology, points)
    halfedgedyads = getHalfEdgeDyads(topology, points)

    r1s, r2s, r3s = getFaceVertices(topology, points)
    facedyads = getFaceDyads(topology, points)
    #### ここまで 118.81 k allocations: 6.495 MiB ####

    for i in eachindex(positions)
        pos = positions[i]

        for j in 1:num_halfedge
            r_i = heads[j]
            r_j = tails[j]
            E_he = halfedgedyads[j]

            pot, acc = getHalfEdgeTerm(r_i, r_j, E_he, pos)
            potentials[i]    += pot
            accelerations[i] += acc
        end

        for j in 1:num_face
            r_i = r1s[j]
            r_j = r2s[j]
            r_k = r3s[j]
            F_f = facedyads[j]

            pot, acc = getFaceTerm(r_i, r_j, r_k, F_f, pos)
            potentials[i]    -= pot
            accelerations[i] -= acc
        end
    end

    G = 6.6740831e-11  # Gravitational constant [m^3/kg/s^2]
    potentials    *= -0.5 * G * ρ
    accelerations *=        G * ρ

    return potentials, accelerations
end

function getHalfEdgeTerm(r_i, r_j, E_he, pos)
    r_he = r_i - pos
    L_he = getEdgeFactor(r_i, r_j, pos)

    acc = E_he * r_he .* L_he  # dyad, vectorの積でallocationが発生する
    pot = r_he ⋅ acc

    return pot, acc
end

function getFaceTerm(r_i, r_j, r_k, F_f, pos)
    r_f = r_i - pos
    ω_f = getFaceFactor(r_i, r_j, r_k, pos)

    acc = F_f * r_f .* ω_f  # vector, dyadの積でallocationが発生する
    pot = r_f ⋅ acc

    return pot, acc
end

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

"""
"""
# function getSurfaceGravity(topology, points, ρ)
#     num_face = nfaces(topology)

#     facecenters = getFaceCenters(topology, points)
#     halfedgedyads = getHalfEdgeDyads(topology, points)
#     facedyads = getFaceDyads(topology, points)

#     potential = zeros(num_face)
#     attraction = zeros(num_face, 3)

#     for i in 1:num_face
#         U, g = getGravity(topology, points, halfedgedyads, facedyads, facecenters[i], ρ)
#         potential[i] = U
#         attraction[i, :] = g
#     end

#     return potential, attraction
# end

"""
    getGravity(topology, points, halfedgedyads, facedyads, pos, ρ=1.) -> U, g

Calculate the gravity field of a constant density polyhedron

# Arguments
- `topology`      : Topology instace of HalfEdges.jl
- `points`        : vertex-list of the polyhedron
- `halfedgedyads` : array of halfedge dyads, E_he
- `facedyads`     : array of face dyads, F_f
- `pos`           : position where you want the gravity
- `ρ`             : constant density of the polyhedron (default: 1.0)
# Returns
- `U` : potential
- `g` : attraction
"""
# function getGravity(topology, points, halfedgedyads, facedyads, pos, ρ)

#     U_e, g_e = sumHalfEdgeTerms(getEdgeVertices(topology, points)..., halfedgedyads, pos)

#     U_f, g_f = sumFaceTerms(getFaceVertices(topology, points)..., facedyads, pos)

#     G = 6.6740831e-11  # Gravitational constant [m^3/kg/s^2]
#     U = -0.5 * G * ρ * (U_e - U_f)
#     g =        G * ρ * (g_e - g_f)

#     return U, g
# end

"""
    sumHalfEdgeTerms(topology, points, halfedgedyads, pos) -> U_he, g_he

Get edge terms of potential and attraction

# Aruguments
- `topology`  : Topology instace of HalfEdges.jl
- `points`        : vertex-list of the polyhedron
- `halfedgedyads` : array of halfedge dyads, E_he
- `pos`           : position where you want the gravity
# Returns
- `U_he` : halfedge term of potential
- `g_he` : halfedge term of attraction
"""
# function sumHalfEdgeTerms(heads, tails, halfedgedyads, pos)
#     num_halfedge = size(halfedgedyads, 1)

#     U_he = 0.
#     g_he = [0., 0., 0.]

#     for i in 1:num_halfedge
#         r_i = heads[i]
#         r_j = tails[i]
#         E_he = halfedgedyads[i]

#         r_he = r_i .- pos
#         L_he = getEdgeFactor(r_i, r_j, pos)

#         U_he += r_he ⋅ (E_he * r_he) * L_he
#         g_he +=        (E_he * r_he) * L_he
#     end
#     return U_he, g_he
# end

"""
    sumFaceTerms(topology, points, facedyads, pos) -> U_f, g_f

Sum face terms over all faces

# Aruguments
- `points`    : vertex-list of the polyhedron
- `faces`     : list of face vertices
- `facedyads` : array of face dyads, F_f
- `pos`       : position where you want the gravity
# Returns
- `U_f` : face term of potential
- `g_f` : face term of attraction
"""
# function sumFaceTerms(r1s, r2s, r3s, facedyads, pos)
#     num_face = size(facedyads, 1)

#     U_f = 0.
#     g_f = [0., 0., 0.]

#     for i in 1:num_face
#         r_i = r1s[i]
#         r_j = r2s[i]
#         r_k = r3s[i]

#         F_f = facedyads[i]

#         r_f = r_i .- pos
#         ω_f = getFaceFactor(r_i, r_j, r_k, pos)

#         U_f += r_f ⋅ (F_f * r_f) * ω_f
#         g_f +=       (F_f * r_f) * ω_f
#     end
#     return U_f, g_f
# end

"""
    getEdgeFactor(r_i, r_j) -> L_e

Calculate a edge factor, L_e

# Arguments
- `r_i` : start-point of the edge
- `r_j` : end-point of the edge
- `pos` : position where you want the gravity (if skipped, the origin is given)
# Returns
- `L_e`: edge factor
"""
function getEdgeFactor(r_i, r_j)
    r_i_norm = norm(r_i)
    r_j_norm = norm(r_j)
    e_ij_norm = norm(r_i .- r_j)

    L_e = r_i_norm + r_j_norm + e_ij_norm
    L_e /= r_i_norm + r_j_norm - e_ij_norm
    L_e = log(L_e)
end

getEdgeFactor(r_i, r_j, pos) = getEdgeFactor(r_i .- pos, r_j .- pos)

"""
    getFaceFactor(r_i, r_j, r_k) -> ω_f

Calculate a face factor, ω_f

# Arguments
- `r_i` : 1st vertex of the face (counter-clockwise)
- `r_j` : 2nd vertex of the face
- `r_k` : 3rd vertex of the face
– `pos` : position where you want the gravity (if skipped, the origin is given)
# Returns
- `ω_f`: face factor
"""
function getFaceFactor(r_i, r_j, r_k)
    r_i_norm = norm(r_i)
    r_j_norm = norm(r_j)
    r_k_norm = norm(r_k)

    y = r_i ⋅ (r_j × r_k)

    x = r_i_norm * r_j_norm * r_k_norm
    x += r_i_norm * (r_j ⋅ r_k)
    x += r_j_norm * (r_k ⋅ r_i)
    x += r_k_norm * (r_i ⋅ r_j)

    ω_f = 2 * atan(y, x)
end

getFaceFactor(r_i, r_j, r_k, pos) = getFaceFactor(r_i .- pos, r_j .- pos, r_k .- pos)

"""
    getHalfEdgeDyads(topology, points) -> halfedgedyads

Get an array of dyads for all halfedges

# Arguments
- `topology` : Topology instace of HalfEdges.jl
- `points`   : vertex-list of the polyhedron
# Return
- halfedgedyads : array of halfedge dyads
"""
function getHalfEdgeDyads(topology, points)
    num_halfedge = nhalfedges(topology)
    halfedgedyads = Array{SMatrix{3,3,Float64}, 1}(undef, num_halfedge)

    for i in 1:num_halfedge
        halfedge = topology.he[i]

        _head = points[halfedge.head]
        _tail = points[HalfEdges.tail(topology, halfedge)]
        edge = normalize(_head - _tail)  # normalized vector from tail to head

        facenormal = HalfEdges.trinormal(topology, points, HalfEdgeHandle(i))
        edgenormal = edge × facenormal
        E_he = getdyad(facenormal, edgenormal)
        halfedgedyads[i] = E_he
    end
    return halfedgedyads
end

"""
    getFaceDyads(topology, points) -> facedyads

Get an array of dyads for all faces

# Arguments
- `topology` : Topology instace of HalfEdges.jl
- `points`   : vertex-list of the polyhedron
# Return
- facedyads : array of face dyads
"""
function getFaceDyads(topology, points)
    num_face = nfaces(topology)
    facedyads = Array{SMatrix{3,3,Float64}, 1}(undef, num_face)
    facenormals = HalfEdges.normals(topology, points)

    for i in 1:num_face
        n̂ = facenormals[i]
        F_f = getdyad(n̂, n̂)
        facedyads[i] = F_f
    end
    return facedyads
end

"""
    getdyad(a, b) -> dyad

Make a dyad from two vectors.
You can simply write: `a * b'`

# Arguments
- `a`, `b`: two vectors
# Return
- `dyad`: dyad
"""
function getdyad(a, b)
    dyad = SA_F64[
        a[1]*b[1] a[1]*b[2] a[1]*b[3];
        a[2]*b[1] a[2]*b[2] a[2]*b[3];
        a[3]*b[1] a[3]*b[2] a[3]*b[3];
    ]
end

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
# This file illustrates the simulation of feedthrough dynamical system in a unity feedback.

using Causal
using Plots

# Construct the model
x0 = ones(1)

@defmodel model begin
    @nodes begin
        gen = FunctionGenerator(readout=sin)
        adder = Adder(signs=(+,-))
        ds = ContinuousLinearSystem(state=x0)
        writer = Writer()
    end
    @branches begin
        gen[1]      =>      adder[1]
        adder       =>      ds
        ds[1]       =>      adder[2]
        ds          =>      writer
    end
end

# Simulate the model
simulate!(model, 0., 0.01, 10.)

# Read simulation data
t, ys = read(getnode(model, :writer).component)

# Compoute simulation error
r = getnode(model,:gen).component.readout.(t)
xa = (x0[1] + 2 / 13) * exp.(-3 / 2 * t) + 3 / 13 * sin.(t) - 2 / 13 * cos.(t)
ya = (xa + r) / 2
er = ys - ya

# Plot results.
p1 = plot(t, ys, label=:simulation)
    plot!(t, ya, label=:analytic)
p2 = plot(t, er, label=:error)
display(plot(p1, p2, layout=(2, 1)))


text:
struct IBPEncoder <: StochasticLayer
    mlp
    ibp
    gauss
end

"""
Forward pass of the encoder that produces the variational distributions.
"""
function (enc::IBPEncoder)(x, d::Int...)
    h = enc.mlp(x)
    dist_kuma, dist_gumbel = enc.ibp(h, d...)
    dist_gauss = (enc.gauss == nothing) ? nothing : enc.gauss(h, d...)
    return dist_kuma, dist_gumbel, dist_gauss
end

function IBPLinearEncoder(k_max::Int, x_dim::Int, h_dim::Int, ibp_enc)
    mlp = Dense(x_dim, h_dim; f=relu)
    return IBPEncoder(mlp, ibp_enc, nothing)
end

function IBPDeepEncoder(k_max::Int, x_dim::Int, h_dim::Int, ibp_enc)
    mlp = Dense(x_dim, h_dim; f=relu)
    gauss = GaussianLogVarLazyDense(h_dim, k_max)
    return IBPEncoder(mlp, ibp_enc, gauss)
end

abstract type AbstractIBPVAE <: NeuralModel end

isdeep(vae::AbstractIBPVAE) = vae.encoder.gauss != nothing

"""
Sample pi in log space. This function also returns nu in log space as the FIRST parameter.
"""
function rand_logpi(dist_nu)
    lognu = logrand(dist_nu)
    logpi = break_logstick_ibp(lognu)
    return lognu, logpi
end

"""
Sample pi in logit space.
"""
function rand_logitpi(dist_nu; return_pi=false)
    logpi = rand_logpi(dist_nu)
    pi = exp.(logpi)
    logitpi = logit.(pi)

    return return_pi ? () : logitpi
end

"""
Make code.

For deep deocder, the code is the Hadamard product of A and Z.
For linear decoder, the code is just Z, as A is treated as the decoder parameter.
"""
function make_code(vae::AbstractIBPVAE, Z, dist_A)
    A = isdeep(vae) ? rand(dist_A) : nothing
    code = isdeep(vae) ? A .* Z : Z
    return code, A
end

"""
This function is used when computing ELBO.
"""
function compute_kl(vae::AbstractIBPVAE, dist_A, dist_nu, dist_Z, logpi, logitpi, logitZ; τ=FT(0.1))
    kl_A = isdeep(vae) ? sum(kldiv(dist_A, BatchNormal(zero(FT), one(FT)))) : 0
    kl_nu = sum(kldiv(dist_nu, BatchBeta(vae.prior.α, one(FT))))
    lq_Z = sum(logpdflogit(dist_Z, logitZ; τ=τ))
    lp_Z = sum(logpdflogit(BatchGumbelBernoulliLogit(logitpi), logitZ; τ=τ))
    kl_Z = lq_Z - lp_Z
    return kl_nu, kl_Z, kl_A
end

"""
This function is used when computing IWAE.
"""
function compute_logratio(vae::AbstractIBPVAE, dist_A, A, dist_nu, nu, pi, dist_Z, Z)
    logratioA = 0
    if isdeep(vae)
        lp_A = sum(logpdf(BatchNormal(zero(FT), one(FT)), A))
        lq_A = sum(logpdf(dist_A, A))
        logratioA = lp_A - lq_A
    end
    lp_nu = sum(logpdf(BatchBeta(vae.prior.α, one(FT)), nu))
    lq_nu = sum(logpdf(dist_nu, nu))
    lp_Z = sum(logpdf(BatchBernoulli(pi), Z))
    lq_Z = sum(logpdf(dist_Z, Z))
    return lp_nu - lq_nu, lp_Z - lq_Z, logratioA
end

"""
Linear Bernoulli decoder for IBP
"""
struct BernoulliIBPLinearDecoder <: StochasticLayer
    A
end

function BernoulliIBPLinearDecoder(k_max, x_dim)
    A = Knet.param(x_dim, k_max; atype=AT{FT,2})
    return BernoulliIBPLinearDecoder(A)
end

function (dec::BernoulliIBPLinearDecoder)(Z)
    k_diff = size(Z, 1) - size(dec.A, 2)
    if k_diff > 0
        dec.A.value = hcat(dec.A.value, Knet.param(size(dec.A, 1), k_diff; atype=AT{FT,2}).value)
    end
    return BatchBernoulliLogit(dec.A[:,1:size(Z,1)] * Z)
end

"""
Linear Gaussian decoder for IBP
"""
struct GaussianIBPLinearDecoder <: StochasticLayer
    sig_X
    A
end

function GaussianIBPLinearDecoder(sig_X, k_max, x_dim)
    A = Knet.param(x_dim, k_max; atype=AT{FT,2})
    return GaussianIBPLinearDecoder(sig_X, A)
end

function (dec::GaussianIBPLinearDecoder)(Z)
    k_diff = size(Z, 1) - size(dec.A, 2)
    if k_diff > 0
        dec.A.value = hcat(dec.A.value, Knet.param(size(dec.A, 1), k_diff; atype=AT{FT,2}).value)
    end
    return BatchNormalLogVar(dec.A[:,1:size(Z,1)] * Z, log(dec.sig_X^2))
end

"""
Write IBP training stats to tbX.
"""
function write_ibp_train_stats!(is_deep, writer, iter, kl_A, kl_nu, kl_Z, lp_x, Z; k_list=nothing, plot_Z=false)
    if writer != nothing
        is_deep && writer."add_scalar"("batch/kl_A", value(kl_A), iter)
        writer."add_scalar"("batch/kl_nu", value(kl_nu), iter)
        writer."add_scalar"("batch/kl_Z", value(kl_Z), iter)
        writer."add_scalar"("batch/nll", -value(lp_x), iter)
        Z_bin = AT == Array ? Z .> 0.5 : Array(Z) .> 0.5
        num_acts = sum(Z_bin; dims=1)
        writer."add_histogram"("batch/num_acts", num_acts, iter)
        num_items = sum(Z_bin; dims=2)
        writer."add_histogram"("batch/num_items", num_items, iter)
        plot_Z && writer."add_image"("batch/Z", Z_bin[:,1:100], iter, dataformats="HW")
        k_list != nothing && writer."add_scalar"("batch/k_avg", mean(k_list), iter)
    end
end


text:
"""
INTERACTIONTrajdata
The trajectory data stored in the original INTERACTION dataset format.
The dataset is a csv file with columns:
    track_id      - Int64 - Representing the id of the agent
    frame_id      - Int64 - Represents the frames in which the agent appears in the video
    timestamp_ms  - Int64 - represents the time the agent appears in the video. The unit is millisecond
    agent_type    - String - Can be person, car, truck and so on
    x             - Float64 - x position, in meter
    y             - Float64 - y position in meter
    vx            - Float64 - x velocity in m/s
    vy            - Float64 - y velocity in m/s
    psi_rad       - Float64 - yaw angle in radian
    length        - Float64 - Length of the vehicle, in meter
    width         - Float64 - Width of the vehicle, in meter

# Example
```julia
tdraw = INTERACTIONTrajdata("vehicle_tracks_000.csv");
```
"""
mutable struct INTERACTIONTrajdata
    df         :: DataFrame
    car2start  :: Dict{Int, Int}         # maps carindex to starting index in the df
    frame2cars :: Dict{Int, Vector{Int}} # maps frame to list of carids in the scene

    function INTERACTIONTrajdata(input_path::String)

        @assert(isfile(input_path))

        df = readtable(input_path, separator=',', header = true)

        car2start = Dict{Int, Int}()
        frame2cars = Dict{Int, Vector{Int}}()

        for (dfind, carid) in enumerate(df[:track_id])
            if !haskey(car2start, carid)
                car2start[carid] = dfind
            end

            frame = convert(Int, df[dfind, :frame_id])
            if !haskey(frame2cars, frame)
                frame2cars[frame] = [carid]
            else
                frame2cars[frame] = push!(frame2cars[frame], carid)
            end
        end

        new(df, car2start, frame2cars)
    end
end

# function: overload nframes from Records
"""
    Records.nframes

# Example
```julia
tdraw = INTERACTIONTrajdata("vehicle_tracks_000.csv")
nframes(tdraw)
```
"""
nframes(trajdata::INTERACTIONTrajdata) = maximum(keys(trajdata.frame2cars))

# function: carsinframe for memory allocation
"""
    carsinframe(trajdata::INTERACTIONTrajdata, frame::Int)
- From the dictionary `frame2cars` which has frame number as key and all vector of car ids in that frame
as value, this function extracts the value corresponding to key given by `frame`

# Example
```julia
# Return array with car ids in frame 1
carsinframe(tdraw,1)
```
"""
carsinframe(trajdata::INTERACTIONTrajdata, frame::Int) = get(trajdata.frame2cars, frame, Int[])

# function: car_df_index - find index of car in trajdata
"""
	function car_df_index

- Given frame and carid, find index of car in trajdata
- Returns 0 if it does not exist (CAVEAT THIS IS COMMENTED OUT)

"""
function car_df_index(trajdata::INTERACTIONTrajdata, carid::Int, frame::Int)
    df = trajdata.df

    lo = trajdata.car2start[carid] # The index in the dataframe where carid first appears
    framestart = df[lo, :frame_id] # The frame number where carid first appears

    retval = 0

    if framestart == frame
        retval = lo
    elseif frame ≥ framestart
        retval = frame - framestart + lo
        #n_frames = df[lo, :n_frames_in_dataset]
        #if retval > lo + n_frames
        #    retval = 0
        #end
    end

    retval
end

"""
    function read_veh_tracks(;roadway)
- Reads the vehicle tracks information from `../dataset/vehicle_tracks_000.csv`
- Modernize to stop using `ListRecord` and `Trajdata`
- Just use Vector of Scenes

# Arguments
- roadway: Used in the vehicle state to create `Frenet` coords for the vehicle on the provided roadway

# Example
```julia
road_ext = make_interaction_roadway_with_extensions()
trajdata = read_veh_tracks(roadway=road_ext)

# Then access a certain framenumber from the trajdata by
scene = trajdata[1]
render([road_ext,scene,IDOverlay(scene=scene)])
```
"""
function read_veh_tracks(;roadway)
    # Read the csv data into tdraw which recasts csv info into a dataframe and also provides two
    # dicts called car2start and frame2cars
    tdraw = INTERACTIONTrajdata(joinpath(@__DIR__,"dataset/vehicle_tracks_000.csv"))
    df = tdraw.df
    vehdefs = Dict{Int, VehicleDef}()

    Trajdata = Vector{Scene}(undef,nframes(tdraw))

    # Initialize the vehicles physically i.e type, length and width
    for (id, dfind) in tdraw.car2start
        # CAUTION: Hardcoding vehicle type to be just car shown by the first argument being 2
        vehdefs[id] = VehicleDef(2, df[dfind, :length], df[dfind, :width])
    end

    # Fill in the vehicle state information in terms of x,y and speed
    for frame in 1 : nframes(tdraw)
        scene = Scene(Entity{VehicleState,VehicleDef,Int64},
            length(carsinframe(tdraw,frame)))
        scene_ind = 0
        for id in carsinframe(tdraw, frame)
            dfind = car_df_index(tdraw, id, frame)

            posG = VecSE2(df[dfind, :x], df[dfind, :y], df[dfind, :psi_rad])
            vx = df[dfind,:vx]
            vy = df[dfind,:vy]
            speed = sqrt(vx*vx + vy*vy)
            push!(scene,Entity(VehicleState(posG, roadway, speed),vehdefs[id],id))
        end

        Trajdata[frame] = scene

    end

    return Trajdata
end

# function: start and end time extraction of vehicles
"""
function extract_timespan(trajscenesvector;offset=1,minlength=20,verbose=false)

- Extract the start time and end time of each vehicle
- Used in sampling a set of vehicles that drive together for a given timespan

# Inputs
- trajscenesvector: Is the output from read_veh_tracks i.e. vector of scenes
- offset from the start

# Returns
- Dict with car id as key and start and end times of the car id

# Examples
```julia
trajscenes = read_veh_tracks(roadway=roadway)
timestamped_trajscenes = extract_timespan(trajscenes)
```
"""
function extract_timespan(trajscenesvector;offset=1,minlength=20,verbose=false)
    index = Dict()
    scene_length = maximum(length(trajscenesvector[i]) for i in 1 : length(trajscenesvector))
    scene = Scene(Entity{VehicleState,VehicleDef,Int64}, scene_length)
    prev, cur = Set(), Set()
    n_frames = length(trajscenesvector)

    # iterate each frame collecting info about the vehicles
    for frame in offset : n_frames - offset
        if verbose > 0
            print("\rframe $(frame) / $(n_frames - offset)")
        end
        cur = Set()
        scene=trajscenesvector[frame]

        # add all the vehicles to the current set
        for veh in scene
            push!(cur, veh.id)
            # insert previously unseen vehicles into the index
            if !in(veh.id, prev)
                index[veh.id] = Dict("ts"=>frame)
            end
        end

        # find vehicles in the previous but not the current frame
        missing = setdiff(prev, cur)
        for id in missing
            # set the final frame for all these vehicles
            index[id]["te"] = frame - 1
        end

        # step forward
        prev = cur
    end

    # at this point, any ids in cur are in the last frame, so add them in
    for id in cur
        index[id]["te"] = n_frames - offset
    end

    # postprocess to remove undesirable trajectories
    for (vehid, infos) in index
        # check for start and end frames
        if !in("ts", keys(infos)) || !in("te", keys(infos))
            if verbose > 0
                println("delete vehid $(vehid) for missing keys")
            end
            delete!(index, vehid)

        # check for start and end frames greater than minlength
        elseif infos["te"] - infos["ts"] < minlength
            if verbose > 0
                println("delete vehid $(vehid) for below minlength")
            end
            delete!(index, vehid)
        end
    end

    return index
end

# Functions: Sampling vehicles that drive together
"""
    function random_sample_from_set_without_replacement(s, n)
Description
    This function samples n values from the set s without replacement, and
    does not work with anything except a set s. Could use statsbase, but want
    to avoid the dependency.
- Used by `sample_multiple_trajdata_vehicle`
Args:
    - s: a set
    - n: number of values to sample
Returns:
    - a subset of the values in s, as a list containing n elements
"""
function random_sample_from_set_without_replacement(s, n)
    @assert length(s) >= n
    sampled = Set()
    for i in 1:n
        cur = rand(s)
        push!(sampled, cur)
        delete!(s, cur)
    end
    return collect(sampled)
end

# Sample vehicles that stay in the scene together
"""
function sample_simultaneous_vehs

- Finds groups of vehicles that stay together on the road for a duration
- Adapted from `ngsim/julia/src/ngsim_utils.jl`

# Inputs
- timestamped_trajscenes: What we get from `extract_timespan` i.e. dict with vehid and ts, te

# Returns
- `egoids`: List of vehicles that drive together with the specified ego vehicle for `offset` num of steps
- `ts`: Start frame number
- `te`: End frame number

# Examples
```julia
a,b,c= sample_simultaneous_vehs(10,50,timestamped_trajscenes,egoid=6)
```
"""
function sample_simultaneous_vehs(
        n_veh::Int,
        offset::Int,
        timestamped_trajscenes;
        max_resamples::Int = 100,
        egoid::Union{Nothing, Int} = nothing,
        verbose::Bool = true,
        rseed::Union{Nothing, Int} = nothing
    )

    if rseed != nothing
        Random.seed!(rseed)
    end

    trajinfos = timestamped_trajscenes

    # if passed in egoid and traj_idx, use those, otherwise, sample
    if egoid == nothing
        # sample the first vehicle and start and end timesteps
        egoid = rand(collect(keys(trajinfos)))
    end

    ts = trajinfos[egoid]["ts"]
    te = trajinfos[egoid]["te"]
    # this sampling assumes ts:te-offset is a valid range
    # this is enforced by the initial computation of the index / trajinfo
    ts = rand(ts:te - offset)
    # after setting the start timestep randomly from the valid range, next
    # update the end timestep to be offset timesteps following it
    # this assume that we just want to simulate for offset timesteps
    te = ts + offset

    # find all other vehicles that have at least 'offset' many steps in common
    # with the first sampled egoid starting from ts. If the number of such
    # vehicles is fewer than n_veh, then resample
    # start with the set containing the first egoid so we don't double count it
    egoids = Set{Int}(egoid)
    for othid in keys(trajinfos)
        oth_ts = trajinfos[othid]["ts"]
        oth_te = trajinfos[othid]["te"]
        # other vehicle must start at or before ts and must end at or after te
        if oth_ts <= ts && te <= oth_te
            push!(egoids, othid)
        end
    end

    # check that there are enough valid ids from which to sample
    if length(egoids) < n_veh
        # if not, resample
        if verbose
            println("WARNING: insuffcient sampling ids in sample_multiple_trajdata_vehicle")
        end
        if max_resamples == 0
            error("ERROR: reached maximum resamples in sample_multiple_trajdata_vehicle")
        else
            return sample_simultaneous_vehs(
                n_veh,
                offset,
                max_resamples=max_resamples - 1,
                verbose=verbose)
        end
    end

    # reaching this point means there are sufficient ids, sample the ones to use
    egoids = random_sample_from_set_without_replacement(egoids, n_veh)

    return egoids, ts, te
end

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
using DTWDT
using Test
using JLD2, Printf

include("../src/functions.jl")
using .DTWDTfunctions

@testset "functions.jl" begin
    #prepare dataset
    maxLag = 80;
    b      = 1;

    fi = @sprintf("%s/testdata.jld2", pwd())
    @eval @load $("$fi"); # data with sine wave shifts of one cycle
    lvec   = (-maxLag:maxLag).*dt; # lag array for plotting below
    npts   = length(u0);            # number of samples
    direction = 1;

    #load true results
    fi = @sprintf("%s/trueresultsfortest.jld2", pwd())
    trueresult = @eval @load $("$fi");
    err_true = eval(trueresult[1])
    dist_true = eval(trueresult[2])
    stbar_true = eval(trueresult[3])
    error_true = eval(trueresult[4])

    err   = computeErrorFunction(u1, u0, npts, maxLag); # cpmpute error function over Lags
    dist  = accumulateErrorFunction(direction, err, npts, maxLag, b); # forward accumulation to make distance function
    stbar = backtrackDistanceFunction( -1*direction, dist, err, -maxLag, b); # find shifts
    error = computeDTWerror(err, stbar, maxLag); #return DTWerror

    @test err == err_true
    @test dist == dist_true
    @test stbar == stbar_true
    @test error == error_true
end


text:
function GRASP_algorithm(beta::Number, prec::Dict{Int, Array}, task_times::Array{Int, 2}, bj::Array{Int,1}, CTS::Constants)
    tasks_by_w = OperationalShippingProblem(beta, task_times, bj, CTS)

    # it_results=Array{Tuple{Number, Number, String, String}, 1}()
    indicators = ["number", "dist"]
    sort_criteria = ["number", "dist"]
    best_makespan = CTS.H
    best_LS = init_ls(CTS)

    for iter = 1:1000
        makespan = CTS.H
        LS = init_ls(CTS)
        QC = init_qc(CTS)
        TIME = init_timer(CTS)
        ind = indicators[rand(1:2)]
        crit = sort_criteria[rand(1:2)]

        it_time = @elapsed while LS.tasks_left > 0 && TIME.period < TIME.horizon_plan
            cranes_status = "idle"
            while cranes_status != "All cranes are busy" && cranes_status != "Next time period" && cranes_status != "LS is completed"
                cranes_status = randomizedConstructionHeuristic(ind, crit, tasks_by_w, prec, bj, LS, TIME, QC, CTS)
            end
            next_time_period(TIME, QC, CTS)
        end

        if check_solution(prec, LS, CTS) == true
            makespan = total_makespan(LS, CTS)
            # push!(it_results, (makespan, it_time, ind, crit))
            if makespan < best_makespan
                best_makespan = makespan
                best_LS = deepcopy(LS)
            end
        end
    end
    # return(it_results, best_makespan, best_LS)
    return(best_makespan, best_LS)
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
using ..SearchProblem
using Random

mutable struct GNP <: SearchProblem.AbstractSearchProblem{Tuple{Int, Int}, Tuple{Int, Int}}
    N_range::AbstractArray{Int}
    density_range::Tuple{Float64,Float64}
    N::Int
    density::Float64
    grid::BitMatrix  #  true = blocked
    rng::MersenneTwister
    function GNP(N_range::AbstractArray{Int}, density_range::Tuple{Real,Real})
        gnp = new(N_range, density_range)
        gnp.rng = MersenneTwister()
        gnp.N = rand(gnp.rng, gnp.N_range)
        gnp.density = rand(gnp.rng, density_range)
        gnp.grid = BitMatrix(undef, gnp.N, gnp.N)
        return gnp
    end
end

function SearchProblem.info(gnp::GNP)
    return Dict(
        :density => gnp.density,
        :grid_side => gnp.N
    )
end

function SearchProblem.obs(gnp::GNP)
    return Float64[gnp.density, gnp.N / maximum(gnp.N_range)]
end

function SearchProblem.reset!(gnp::GNP)
    gnp.N = rand(gnp.rng, gnp.N_range)
    density = gnp.density_range[1] + (gnp.density_range[2] - gnp.density_range[1]) * rand(gnp.rng)
    I = gnp.N ÷ 20
    gnp.grid = rand(gnp.rng, gnp.N, gnp.N) .< (density / I)
    for i in 1:I
        gnp.grid = gnp.grid .| gnp.grid[:, vcat(2:gnp.N, 1)]
    end
    gnp.density = sum(gnp.grid) / (gnp.N^2)
end

@inline function SearchProblem.start_state(::GNP)::Tuple{Int, Int}
    return (1, 1)
end

function SearchProblem.cost(::GNP, ::Tuple{Int, Int}, ::Tuple{Int, Int}, ::Tuple{Int, Int})::Float64
    return 1.0
end

@inline function SearchProblem.goal_test(gnp::GNP, state::Tuple{Int, Int})::Bool
    return state == (gnp.N, gnp.N)
end

function SearchProblem.heuristic(gnp::GNP, state::Tuple{Int, Int})::Float64
    return (state .- gnp.N) .|> abs |> sum
end

function SearchProblem.successors(gnp::GNP, state::Tuple{Int, Int})
    @inbounds cur_city = state[end]
    available_actions = Iterators.filter(((0, 1), (0, -1), (1, 0), (-1, 0))) do a
        return all(1 .<= state .+ a  .<= gnp.N)  &&  !gnp.grid[(state .+ a)...]
    end
    return Iterators.map(available_actions) do action
        next_state = state .+ action
        return next_state, action
    end
end


text:
using Documenter
using ParticipationData

makedocs(
    modules = [ParticipationData],
    sitename = "ParticipationData.jl",
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(repo = "github.com/laschuet/ParticipationData.jl.git")


text:
using Pkg.Artifacts
using Flux
using BSON

"""
    cats_dogs_data()

A convenience function to load the cats and dogs toy dataset.
"""
# Cats and dogs:
function cats_dogs_data()
    data_dir = artifact"cats_dogs_data"
    data = BSON.load(joinpath(data_dir,"cats_dogs_data.bson"),@__MODULE__)[:data]
    X, y = (data[:X], data[:y])
    return X, y
end

"""
    cats_dogs_model()

A convenience function to load the pre-trained MLP to classify cats and dogs.
"""
function cats_dogs_model()
    data_dir = artifact"cats_dogs_model"
    model = BSON.load(joinpath(data_dir,"cats_dogs_model.bson"),@__MODULE__)[:model]
    return model
end

"""
    cats_dogs_laplace()

A convenience function to load the pre-trained MLP with Laplace approximation to classify cats and dogs.
"""
function cats_dogs_laplace()
    data_dir = artifact"cats_dogs_laplace"
    la = BSON.load(joinpath(data_dir,"cats_dogs_laplace.bson"),@__MODULE__)[:la]
    return la
end

# MNIST:
"""
    mnist_data()

A convenience function to load MNIST training data.
"""
function mnist_data()
    data_dir = artifact"mnist_data"
    data = BSON.load(joinpath(data_dir,"mnist_data.bson"),@__MODULE__)[:data]
    X, ys = (data[:X], data[:ys])
    return X, ys
end

"""
    mnist_model()

A convenience function to load the pre-trained MLP for MNIST training data.
"""
function mnist_model()
    data_dir = artifact"mnist_model"
    model = BSON.load(joinpath(data_dir,"mnist_model.bson"),@__MODULE__)[:model]
    return testmode!(model)
end

"""
    mnist_ensemble()

A convenience function to load the pre-trained deep ensemble of MLPs for MNIST training data.
"""
function mnist_ensemble()
    data_dir = joinpath(artifact"mnist_ensemble","mnist_ensemble")
    model_files = Base.Filesystem.readdir(data_dir)
    ensemble = []
    for file in model_files
        model = BSON.load(joinpath(data_dir,file),@__MODULE__)[:model]
        ensemble = vcat(ensemble, testmode!(model))
    end
    return ensemble
end

using Random
"""
    toy_data_linear(N=100)

A convenience function to load linearly separable synthetic data.

# Examples

```julia-repl
toy_data_linear()
```

"""
function toy_data_linear(N=100)
    # Number of points to generate.
    M = round(Int, N / 2)
    Random.seed!(1234)

    # Generate artificial data.
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt1s = Array([[x1s[i] + 0.5; x2s[i] + 0.5] for i = 1:M])
    xt0s = Array([[x1s[i] - 5; x2s[i] - 5] for i = 1:M])

    # Store all the data for later.
    xs = [xt1s; xt0s]
    ts = [ones(M); zeros(M)];
    return xs, ts
end

using Random
"""
    toy_data_non_linear(N=100)

A convenience function to load synthetic data that are not linearly separable.

# Examples

```julia-repl
toy_data_non_linear()
```

"""
function toy_data_non_linear(N=100)
    # Number of points to generate.
    M = round(Int, N / 4)
    Random.seed!(1234)

    # Generate artificial data.
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt1s = Array([[x1s[i] + 0.5; x2s[i] + 0.5] for i = 1:M])
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    append!(xt1s, Array([[x1s[i] - 5; x2s[i] - 5] for i = 1:M]))

    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt0s = Array([[x1s[i] + 0.5; x2s[i] - 5] for i = 1:M])
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    append!(xt0s, Array([[x1s[i] - 5; x2s[i] + 0.5] for i = 1:M]))

    # Store all the data for later.
    xs = [xt1s; xt0s]
    ts = [ones(2*M); zeros(2*M)];
    return xs, ts
end

using Random
"""
    toy_data_multi(N=100)

A convenience function to load multi-class synthetic data.

# Examples

```julia-repl
toy_data_multi()
```

"""
function toy_data_multi(N=100)
    # Number of points to generate.
    M = round(Int, N / 4)
    Random.seed!(1234)

    # Generate artificial data.
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt1s = Array([[x1s[i] + 1; x2s[i] + 1] for i = 1:M])
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    append!(xt1s, Array([[x1s[i] - 7; x2s[i] - 7] for i = 1:M]))

    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt0s = Array([[x1s[i] + 1; x2s[i] - 7] for i = 1:M])
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    append!(xt0s, Array([[x1s[i] - 7; x2s[i] + 1] for i = 1:M]))

    # Store all the data for later.
    xs = [xt1s; xt0s]
    ts = [ones(M); ones(M).*2; ones(M).*3; ones(M).*4];
    return xs, ts
end


text:
# Copyright 2021 Gravifer
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test

include(normpath(joinpath(@__DIR__, "../src/CFD2021Projects.jl")))
@testset "misc_utils" begin
    include(normpath(joinpath(@__DIR__, "../src/misc-util.jl")))

    @testset "tuplejoin" begin
        using .__CFD2021__misc_util__: tuplejoin
        @test tuplejoin((1, 2)) == (1, 2)
        temp_v = ([0,0.5,1], [0,0.5,1], [0,0.5,1], [0,0.5,1])
        temp_t = (axes(v) for v in temp_v)
        @test tuplejoin(temp_t...) == (Base.OneTo(3), Base.OneTo(3), Base.OneTo(3), Base.OneTo(3))
    end # @testset "tuplejoin"
end # @testset "misc_utils"

@testset "Project-1" begin
    include(normpath(joinpath(@__DIR__, "../Project-1/test/runtests.jl")))
end # @testset "Project-1"


text:
using BinaryBuilder, ObjectFile

verbose = "--verbose" in ARGS

# GR already contains a bunch of tarballs for us, so we basically just download
# and rename them for our own nefarious purposes.

name = "GR"
version = v"0.41.0"
name_mapping = Dict(
    Linux(:x86_64, libc=:glibc) => ("Linux-x86_64", "516bbd70640f4f49df5968d9d85513c6f5fa923681f59c1876d54c3b055bf2c3"),
    Linux(:i686, libc=:glibc) => ("Linux-i386", "095b98bef33f3e140c7f8ed0008685b9c79925ee5c45545242e51fce837dd7a9"),
    Windows(:x86_64) => ("Windows-x86_64", "ac7437bcb067cefa17c54d6851d18d4f8220dee8bb4177e4f7c5761965c0d6cc"),
    Windows(:i686) => ("Windows-i686", "5e9ef15f0a53746fa6db965959458992314f59d33c013f693b8be3e093420e4e"),
    MacOS(:x86_64) => ("Darwin-x86_64", "1aaa4c6cfef141a3298518178e8f8dc0f73b1ac8ab5ef826c59158ab6a590500"),
)

# Downoad, unpack, extract, then repackage each of these guys
mkpath(joinpath(@__DIR__, "build"))
mkpath(joinpath(@__DIR__, "products"))
product_hashes = Dict()
for (platform, (suffix, hash)) in name_mapping
    extract_dir = joinpath(@__DIR__, "build", triplet(platform))
    rm(extract_dir; force=true, recursive=true)
    tarball_path = joinpath(@__DIR__, "build", "$(name)-v$(version)-$(triplet(platform)).tar.gz")
    BinaryBuilder.download_verify_unpack("https://github.com/sciapp/gr/releases/download/v$(version)/gr-$(version)-$(suffix).tar.gz", hash, extract_dir; tarball_path=tarball_path, ignore_existence=true, force=true, verbose=verbose)

    if platform isa MacOS
        # If we're dealing with a MacOS build, rename the top-level `.so` objects to `.dylib`
        libdir = joinpath(extract_dir, "gr", "lib")
        symlink("libGR.so", joinpath(libdir, "libGR.dylib"))
        symlink("libGR3.so", joinpath(libdir, "libGR3.dylib"))
    end

    tarball_path, hash = BinaryBuilder.package(Prefix(joinpath(extract_dir, "gr")), joinpath(@__DIR__, "products", name), version; platform=platform, verbose=verbose, force=true)
    product_hashes[triplet(platform)] = (basename(tarball_path), hash)
end

products = [
    LibraryProduct(Prefix(pwd()), "libGR", :libGR),
    LibraryProduct(Prefix(pwd()), "libGR3", :libGR3),
]

bin_path = "https://github.com/$(BinaryBuilder.get_repo_name())/releases/download/$(BinaryBuilder.get_tag_name())"
BinaryBuilder.print_buildjl(@__DIR__, name, version, products, product_hashes, bin_path)


text:
#=
This file establishes the agent-space interaction API.
=#
export move_agent!, add_agent!, add_agent_single!, add_agent_pos!,
move_agent_single!, kill_agent!, genocide!, nextid

#######################################################################################
# Killing agents
#######################################################################################
"""
    kill_agent!(agent::AbstractAgent, model::ABM)

Remove an agent from model, and from the space if the model has a space.
"""
function kill_agent!(agent::AbstractAgent, model::ABM{A, S}) where {A, S<:AbstractSpace}
  agentnode = coord2vertex(agent.pos, model)
   # remove from the space
  splice!(agent_positions(model)[agentnode],
          findfirst(a->a==agent.id, agent_positions(model)[agentnode]))
  delete!(model.agents, agent.id)
  return model
end

function kill_agent!(agent::A, model::ABM{A, Nothing}) where A
  delete!(model.agents, agent.id)
end

kill_agent!(id::Integer, model) = kill_agent!(model[id], model)

"""
    genocide!(model::ABM)
Kill all the agents of the model.
"""
genocide!(model::ABM) = for (i, a) in model.agents; kill_agent!(a, model); end

"""
    genocide!(model::ABM, n::Int)
Kill the agents of the model whose IDs are larger than n.
"""
function genocide!(model::ABM, n::Int)
    for (k, v) in model.agents
        k > n && kill_agent!(v, model)
    end
end

"""
    genocide!(model::ABM, f::Function)
Kill all agents where the function `f(agent)` returns `true`.
"""
function genocide!(model::ABM, f::Function)
    for (k, v) in model.agents
        f(v) && kill_agent!(v, model)
    end
end

#######################################################################################
# Moving agents
#######################################################################################
"""
    move_agent!(agent::A [, pos], model::ABM{A, <: DiscreteSpace}) → agent

Add `agentID` to the new position `pos` (or a random one if `pos` is not given)
in the model and remove it from the old position
(also update the agent to have the new position).
`pos` must be the appropriate position type depending on the space type.
"""
function move_agent!(agent::AbstractAgent, pos::Tuple, model::ABM)
  nodenumber = coord2vertex(pos, model)
  move_agent!(agent, nodenumber, model)
end

function move_agent!(agent::AbstractAgent, pos::Integer, model::ABM)
  # remove agent from old position
  if typeof(agent.pos) <: Tuple
    oldnode = coord2vertex(agent.pos, model)
    splice!(model.space.agent_positions[oldnode], findfirst(a->a==agent.id, model.space.agent_positions[oldnode]))
    agent.pos = vertex2coord(pos, model)  # update agent position
  else
    splice!(model.space.agent_positions[agent.pos], findfirst(a->a==agent.id, model.space.agent_positions[agent.pos]))
    agent.pos = pos
  end
  push!(model.space.agent_positions[pos], agent.id)
  return agent
end

function move_agent!(agent::AbstractAgent, model::ABM)
  nodenumber = rand(1:nv(model.space))
  move_agent!(agent, nodenumber, model)
  return agent
end

"""
    move_agent_single!(agent::AbstractAgent, model::ABM) → agent

Move agent to a random node while respecting a maximum of one agent
per node. If there are no empty nodes, the agent wont move.
Only valid for non-continuous spaces.
"""
function move_agent_single!(agent::AbstractAgent, model::ABM)
  # TODO: this inefficient
  empty_cells = [i for i in 1:length(model.space.agent_positions) if length(model.space.agent_positions[i]) == 0]
  if length(empty_cells) > 0
    random_node = rand(empty_cells)
    move_agent!(agent, random_node, model)
  end
  return agent
end

#######################################################################################
# Adding agents
#######################################################################################
# TODO: the Source code of `add_agent!` for discrete space is overly complex.
# Similarly with the continuous
# space case, `add_agent_pos!` should be defined first and be the low-level method
# all other methods call. The other `add_agent!` methods simply prepare an agent with
# appropriate position and give it to `add_agent_pos!`

"""
    add_agent_pos!(agent::AbstractAgent, model::ABM) → agent
Add the agent to the `model` at the agent's own position.
"""
function add_agent_pos!(agent::AbstractAgent, model::ABM)
  if :pos ∈ fieldnames(typeof(agent))
    return add_agent!(agent, agent.pos, model)
  else
    return add_agent!(agent, model)
  end
end

"""
    add_agent!(agent::AbstractAgent [, position], model::ABM) → agent

Add the `agent` to the `position` in the space and to the list of agents.
If `position` is not given, the `agent` is added to a random position.
The `agent`'s position is always updated to match `position`, and therefore for `add_agent!`
the position of the `agent` is meaningless. Use [`add_agent_pos!`](@ref) to use
the `agent`'s position.
"""
function add_agent!(agent::A, pos::Tuple, model::ABM{A, <: DiscreteSpace}) where {A}
  # node number from x, y, z coordinates
  nodenumber = coord2vertex(pos, model)
  add_agent!(agent, nodenumber, model)
end

function add_agent!(agent::A, pos::Integer, model::ABM{A, <: DiscreteSpace}) where {A}
  push!(model.space.agent_positions[pos], agent.id)
  model[agent.id] = agent
  # update agent position
  agent.pos = correct_pos_type(pos, model)
  return agent
end

function add_agent!(agent::A, model::ABM{A, <: DiscreteSpace}) where {A}
  nodenumber = rand(1:nv(model.space))
  add_agent!(agent, nodenumber, model)
end

function add_agent!(agent::A, model::ABM{A, <: Nothing}) where {A}
  model[agent.id] = agent
  return agent
end

function add_agents!(agents::Vector{A}, args...) where {A}
  for agent in agents
    add_agent!(agent, args...)
  end
end

"""
    add_agent!([pos,] model::ABM, args...; kwargs...)
Create and add a new agent to the model by constructing an agent of the
type of the `model`. Propagate all *extra* positional arguments and
keyword arguemts to the agent constructor.

Notice that this function takes care of setting the agent's id *and* position and thus
`args...` and `kwargs...` are propagated to other fields the agent has.

Optionally provide a position to add the agent to as first argument.

## Example
```julia
using Agents
mutable struct Agent <: AbstractAgent
    id::Int
    pos::Int
    w::Float64
    k::Bool
end
Agent(id; w, k) = Agent(id, w, k) # keyword constructor
model = ABM(Agent, GraphSpace(complete_digraph(5)))

add_agent!(model, 1, 0.5, true) # incorrect: id is set internally
add_agent!(model, 0.5, true) # correct: weight becomes 0.5
add_agent!(5, model, 0.5, true) # add at node 5
add_agent!(model; w = 0.5, k =true) # use keywords: weight becomes 0.5
```
"""
function add_agent!(node, model::ABM{A, <: DiscreteSpace}, properties...; kwargs...) where {A}
    id = nextid(model)
    cnode = correct_pos_type(node, model)
    agent = A(id, cnode, properties...; kwargs...)
    add_agent!(agent, cnode, model)
end

function add_agent!(model::ABM{A, Nothing}, properties...; kwargs...) where {A}
  @assert model.space == nothing
  id = nextid(model)
  model[id] = A(id, properties...; kwargs...)
  return model[id]
end

function add_agent!(model::ABM{A, Nothing}, ::Type{B}, properties...; kwargs...) where {A, B <: A}
  @assert model.space == nothing
  id = nextid(model)
  model[id] = B(id, properties...; kwargs...)
  return model[id]
end

function add_agent!(model::ABM{A, S}, properties...; kwargs...) where {A, S<:AbstractSpace}
  id = nextid(model)
  n = rand(1:nv(model))
  cnode = correct_pos_type(n, model)
  model[id] = A(id, cnode, properties...; kwargs...)
  push!(model.space.agent_positions[n], id)
  return model[id]
end

"""
    nextid(model::ABM) → id
Return a valid `id` for creating a new agent with it.
"""
nextid(model::ABM) = isempty(model.agents) ? 1 : maximum(keys(model.agents)) + 1

"""
    add_agent_single!(agent::A, model::ABM{A, <: DiscreteSpace}) → agent

Add agent to a random node in the space while respecting a maximum one agent per node.
This function throws a warning if no empty nodes remain.
"""
function add_agent_single!(agent::A, model::ABM{A, <: DiscreteSpace}) where {A}
  msa = model.space.agent_positions
  empty_cells = [i for i in 1:length(msa) if length(msa[i]) == 0]
  if length(empty_cells) > 0
    random_node = rand(empty_cells)
    add_agent!(agent, random_node, model)
  else
    @warn "No empty nodes found for `add_agent_single!`."
  end
  return agent
end

"""
    add_agent_single!(model::ABM{A, <: DiscreteSpace}, properties...; kwargs...)
Same as `add_agent!(model, properties...)` but ensures that it adds an agent
into a node with no other agents (does nothing if no such node exists).
"""
function add_agent_single!(model::ABM{A, <: DiscreteSpace}, properties...; kwargs...) where {A}
  msa = model.space.agent_positions
  id = nextid(model)
  empty_cells = [i for i in 1:length(msa) if length(msa[i]) == 0]
  if length(empty_cells) > 0
    random_node = rand(empty_cells)
    cnode = correct_pos_type(node, model)
    agent = A(id, cnode, properties...; kwargs...)
    add_agent!(agent, random_node, model)
    return agent
  end
end


text:
mutable struct CounterfactualExplanation
    x::AbstractArray
    target::Number
    target_encoded::Union{Number, AbstractVector, Nothing}
    x′::AbstractArray
    data::DataPreprocessing.CounterfactualData
    M::Models.AbstractFittedModel
    generator::Generators.AbstractGenerator
    params::Dict
    search::Union{Dict,Nothing}
end

using .Counterfactuals
function Base.show(io::IO, z::CounterfactualExplanation)

    printstyled(io, "Factual: ", bold=true)
    println(io, "x=$(z.x), y=$(Counterfactuals.factual_label(z)), p=$(Counterfactuals.factual_probability(z))")
    printstyled(io, "Target: ", bold=true)
    println(io, "target=$(z.target), γ=$(z.params[:γ])")

    if !isnothing(z.search)
        printstyled(io, "Counterfactual outcome: ", bold=true)
        println(io, "x′=$(z.x′), y′=$(Counterfactuals.counterfactual_label(z)), p′=$(Counterfactuals.counterfactual_probability(z))")
        printstyled(io, "Converged: $(Counterfactuals.converged(z) ? "✅"  : "❌") ", bold=true)
        println("after $(Counterfactuals.total_steps(z)) steps.")
    else
        @info "Search not yet initatiated."
    end

end

"""
    CounterfactualExplanation(
        x::Union{AbstractArray,Int},
        target::Union{AbstractFloat,Int},
        data::CounterfactualData,
        M::Models.AbstractFittedModel,
        generator::Generators.AbstractGenerator,
        γ::AbstractFloat,
        T::Int
    )

Outer method to construct a `CounterfactualExplanation` structure.
"""
# Outer constructor method:
function CounterfactualExplanation(
    x::Union{AbstractArray,Int},
    target::Union{AbstractFloat,Int},
    data::CounterfactualData,
    M::Models.AbstractFittedModel,
    generator::Generators.AbstractGenerator,
    γ::AbstractFloat,
    T::Int
)
    # Factual:
    x = typeof(x) == Int ? select_factual(data, x) : x
    # Counterfactual:
    x′ = copy(x)  # start from factual

    # Parameters:
    params = Dict(
        :γ => γ,
        :T => T,
        :mutability => DataPreprocessing.mutability_constraints(data)
    )

    return CounterfactualExplanation(x, target, nothing, x′, data, M, generator, params, nothing)

end

# Convenience methods:

# 0) Utils
"""
    output_dim(counterfactual_explanation::CounterfactualExplanation)

A convenience method that computes the output dimension of the predictive model.
"""
output_dim(counterfactual_explanation::CounterfactualExplanation) = size(Models.probs(counterfactual_explanation.M, counterfactual_explanation.x))[1]

using Flux
"""
    encode_target(counterfactual_explanation::CounterfactualExplanation)

A convenience method to encode the target variable, if necessary.
"""
function encode_target(counterfactual_explanation::CounterfactualExplanation)
    out_dim = output_dim(counterfactual_explanation)
    target = counterfactual_explanation.target
    return out_dim > 1 ? Flux.onehot(target, 1:out_dim) : target
end

# 1) Factual values
"""
    factual(counterfactual_explanation::CounterfactualExplanation)

A convenience method to get the factual value.
"""
factual(counterfactual_explanation::CounterfactualExplanation) = counterfactual_explanation.x

"""
    factual_probability(counterfactual_explanation::CounterfactualExplanation)

A convenience method to compute the class probabilities of the factual value.
"""
factual_probability(counterfactual_explanation::CounterfactualExplanation) = Models.probs(counterfactual_explanation.M, counterfactual_explanation.x)

"""
    factual_label(counterfactual_explanation::CounterfactualExplanation)

A convenience method to get the predicted label associated with the factual value.
"""
function factual_label(counterfactual_explanation::CounterfactualExplanation)
    p = factual_probability(counterfactual_explanation)
    out_dim = size(p)[1]
    y = out_dim == 1 ? round(p[1]) : Flux.onecold(p,1:out_dim)
    return y
end

# 2) Counterfactual values:
"""
    initialize!(counterfactual_explanation::CounterfactualExplanation)

A subroutine that intializes the counterfactual search.
"""
function initialize!(counterfactual_explanation::CounterfactualExplanation)

    # Encode target:
    counterfactual_explanation.target_encoded = encode_target(counterfactual_explanation)

    # Initialize search:
    counterfactual_explanation.search = Dict(
        :iteration_count => 1,
        :times_changed_features => zeros(length(counterfactual_explanation.x)),
        :path => [counterfactual_explanation.x′],
        :terminated => threshold_reached(counterfactual_explanation),
        :converged => threshold_reached(counterfactual_explanation),
        :other => nothing
    )

    if counterfactual_explanation.search[:terminated]
        @info "Factual already in target class and probability exceeds threshold γ."
    end

end

"""
    counterfactual(counterfactual_explanation::CounterfactualExplanation)

A convenience method to get the counterfactual value.
"""
counterfactual(counterfactual_explanation::CounterfactualExplanation) = counterfactual_explanation.x′

"""
    counterfactual_probability(counterfactual_explanation::CounterfactualExplanation)

A convenience method to compute the class probabilities of the counterfactual value.
"""
counterfactual_probability(counterfactual_explanation::CounterfactualExplanation) = Models.probs(counterfactual_explanation.M, counterfactual_explanation.x′)

"""
    counterfactual_label(counterfactual_explanation::CounterfactualExplanation)

A convenience method to get the predicted label associated with the counterfactual value.
"""
function counterfactual_label(counterfactual_explanation::CounterfactualExplanation)
    p = counterfactual_probability(counterfactual_explanation)
    out_dim = size(p)[1]
    y = out_dim == 1 ? round(p[1]) : Flux.onecold(p,1:out_dim)
    return y
end

# 3) Search related methods:
"""
    terminated(counterfactual_explanation::CounterfactualExplanation)

A convenience method to determine if the counterfactual search has terminated.
"""
terminated(counterfactual_explanation::CounterfactualExplanation) = counterfactual_explanation.search[:terminated]

"""
    converged(counterfactual_explanation::CounterfactualExplanation)

A convenience method to determine if the counterfactual search has converged.
"""
converged(counterfactual_explanation::CounterfactualExplanation) = counterfactual_explanation.search[:converged]

"""
    total_steps(counterfactual_explanation::CounterfactualExplanation)

A convenience method that returns the total number of steps of the counterfactual search.
"""
total_steps(counterfactual_explanation::CounterfactualExplanation) = counterfactual_explanation.search[:iteration_count]

"""
    path(counterfactual_explanation::CounterfactualExplanation)

A convenience method that returns the entire counterfactual path.
"""
path(counterfactual_explanation::CounterfactualExplanation) = counterfactual_explanation.search[:path]

"""
    target_probs(counterfactual_explanation::CounterfactualExplanation, x::Union{AbstractArray, Nothing}=nothing)

Returns the predicted probability of the target class for `x`. If `x` is `nothing`, the predicted probability corresponding to the counterfactual value is returned.
"""
function target_probs(counterfactual_explanation::CounterfactualExplanation, x::Union{AbstractArray, Nothing}=nothing)

    p = !isnothing(x) ? Models.probs(counterfactual_explanation.M, x) : counterfactual_probability(counterfactual_explanation)
    target = counterfactual_explanation.target

    if length(p) == 1
        h(x) = ifelse(x==-1,0,x)
        if target ∉ [0,1] && target ∉ [-1,1]
            throw(DomainError("For binary classification expecting target to be in {0,1} or {-1,1}."))
        end
        # If target is binary (i.e. outcome 1D from sigmoid), compute p(y=0):
        p = vcat(1.0 .- p, p)
        # Choose first (target+1) row if target=0, second row (target+1) if target=1:
        p_target = p[Int(h(target)+1),:]
    else
        if target < 1 || target % 1 !=0
            throw(DomainError("For multi-class classification expecting `target` ∈ ℕ⁺, i.e. {1,2,3,...}."))
        end
        # If target is multi-class, choose corresponding row (e.g. target=2 -> row 2)
        p_target = p[Int(target),:]
    end
    return p_target
end

"""
    apply_mutability(Δx′::AbstractArray, counterfactual_explanation::CounterfactualExplanation)

A subroutine that applies mutability constraints to the proposed vector of feature perturbations.
"""
function apply_mutability(Δx′::AbstractArray, counterfactual_explanation::CounterfactualExplanation)

    mutability = counterfactual_explanation.params[:mutability]
    # Helper functions:
    both(x) = x
    increase(x) = ifelse(x<0,0,x)
    decrease(x) = ifelse(x>0,0,x)
    none(x) = 0
    cases = (both = both, increase = increase, decrease = decrease, none = none)

    # Apply:
    Δx′ = [getfield(cases, mutability[d])(Δx′[d]) for d in 1:length(Δx′)]

    return Δx′

end

"""
    threshold_reached(counterfactual_explanation::CounterfactualExplanation)

A convenience method that determines of the predefined threshold for the target class probability has been reached.
"""
threshold_reached(counterfactual_explanation::CounterfactualExplanation) = target_probs(counterfactual_explanation)[1] >= counterfactual_explanation.params[:γ]

"""
    steps_exhausted(counterfactual_explanation::CounterfactualExplanation)

A convenience method that checks if the number of maximum iterations has been exhausted.
"""
steps_exhausted(counterfactual_explanation::CounterfactualExplanation) = counterfactual_explanation.search[:iteration_count] == counterfactual_explanation.params[:T]

"""
    get_counterfactual_state(counterfactual_explanation::CounterfactualExplanation)

A subroutine that is used to take a snapshot of the current counterfactual search state. This snapshot is passed to the counterfactual generator.
"""
function get_counterfactual_state(counterfactual_explanation::CounterfactualExplanation)
    counterfactual_state = Generators.CounterfactualState(
        counterfactual_explanation.x,
        counterfactual_explanation.target_encoded,
        counterfactual_explanation.x′,
        counterfactual_explanation.M,
        counterfactual_explanation.params,
        counterfactual_explanation.search
    )
    return counterfactual_state
end

"""
    update!(counterfactual_explanation::CounterfactualExplanation)

An important subroutine that updates the counterfactual explanation. It takes a snapshot of the current counterfactual search state and passes it to the generator. Based on the current state the generator generates perturbations. Various constraints are then applied to the proposed vector of feature perturbations. Finally, the counterfactual search state is updated.
"""
function update!(counterfactual_explanation::CounterfactualExplanation)

    counterfactual_state = get_counterfactual_state(counterfactual_explanation)

    # Generate peturbations:
    Δx′ = Generators.generate_perturbations(counterfactual_explanation.generator, counterfactual_state)
    Δx′ = apply_mutability(Δx′, counterfactual_explanation)
    Δx′ = reshape(Δx′, size(counterfactual_explanation.x′))
    x′ = counterfactual_explanation.x′ + Δx′
    x′ = DataPreprocessing.apply_domain_constraints(counterfactual_explanation.data, x′)
    counterfactual_explanation.x′ = x′ # update counterfactual

    # Updates:
    counterfactual_explanation.search[:path] = [counterfactual_explanation.search[:path]..., counterfactual_explanation.x′]
    counterfactual_explanation.search[:mutability] = Generators.mutability_constraints(counterfactual_explanation.generator, counterfactual_state)
    counterfactual_explanation.search[:times_changed_features] += reshape(Δx′ .!= 0, size(counterfactual_explanation.search[:times_changed_features])) # update number of times feature has been changed
    counterfactual_explanation.search[:iteration_count] += 1 # update iteration counter
    counterfactual_explanation.search[:converged] = threshold_reached(counterfactual_explanation)
    counterfactual_explanation.search[:terminated] = counterfactual_explanation.search[:converged] || steps_exhausted(counterfactual_explanation) || Generators.conditions_satisified(counterfactual_explanation.generator, counterfactual_state)
end


text:
# Autogenerated wrapper script for MPFR_jll for x86_64-w64-mingw32
export libmpfr

using GMP_jll
JLLWrappers.@generate_wrapper_header("MPFR")
JLLWrappers.@declare_library_product(libmpfr, "libmpfr-6.dll")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libmpfr,
        "bin\\libmpfr-6.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Dates: Period, DatePeriod, TimePeriod, floorceil

function Base.floor(lzdt::LaxZonedDateTime, p::DatePeriod)
    return LaxZonedDateTime(floor(DateTime(lzdt), p), timezone(lzdt))
end

function Base.floor(lzdt::LaxZonedDateTime, p::TimePeriod)
    # Rounding non-representable dates doesn't work.
    !isrepresentable(lzdt) && (return LaxZonedDateTime())

    local_dt = DateTime(lzdt)
    local_dt_floored = floor(local_dt, p)

    if local_dt == local_dt_floored
        lzdt
    elseif isvalid(lzdt)
        # Rounding is done using the current fixed offset to avoid transitional ambiguities.
        utc_dt_floored = local_dt_floored - lzdt.zone.offset
        LaxZonedDateTime(ZonedDateTime(utc_dt_floored, timezone(lzdt); from_utc=true))
    else
        LaxZonedDateTime()
    end
end

function Base.ceil(lzdt::LaxZonedDateTime, p::DatePeriod)
    return LaxZonedDateTime(ceil(DateTime(lzdt), p), timezone(lzdt))
end

# TODO: Additional performance gains can be made for round

function Base.round(lzdt::LaxZonedDateTime, p::DatePeriod, r::RoundingMode{:NearestTiesUp})
    f, c = floorceil(lzdt, p)

    if f == c
        f
    elseif isvalid(f) && isvalid(c)
        local_dt = DateTime(lzdt)
        local_dt - DateTime(f) < DateTime(c) - local_dt ? f : c
    else
        LaxZonedDateTime()
    end
end

function Base.round(lzdt::LaxZonedDateTime, p::TimePeriod, r::RoundingMode{:NearestTiesUp})
    f, c = floorceil(lzdt, p)

    if f == c
        f
    elseif isvalid(f) && isvalid(c)
        utc_dt = DateTime(lzdt, UTC)
        utc_dt - DateTime(f, UTC) < DateTime(c, UTC) - utc_dt ? f : c
    else
        LaxZonedDateTime()
    end
end

"""
    floor(lzdt::LaxZonedDateTime, p::Period) -> LaxZonedDateTime
    floor(lzdt::LaxZonedDateTime, p::Type{Period}) -> LaxZonedDateTime

Returns the nearest `LaxZonedDateTime` less than or equal to `zdt` at resolution `p`. The
result will be in the same time zone as `lzdt`.

For convenience, `p` may be a type instead of a value: `floor(lzdt, Dates.Hour)` is a
shortcut for `floor(lzdt, Dates.Hour(1))`.

`VariableTimeZone` transitions are handled as for `round`.

### Examples

The `America/Winnipeg` time zone transitioned from Central Standard Time (UTC-6:00) to
Central Daylight Time (UTC-5:00) on 2016-03-13, moving directly from 01:59:59 to 03:00:00.

```julia
julia> lzdt = LaxZonedDateTime(DateTime(2016, 3, 13, 1, 45), TimeZone("America/Winnipeg"))
2016-03-13T01:45:00-06:00

julia> floor(lzdt, Dates.Hour)
2016-03-13T01:00:00-06:00

julia> floor(lzdt, Dates.Day)
2016-03-13T00:00:00-06:00
```

The `Asia/Colombo` time zone revised the definition of Lanka Time from UTC+6:30 to UTC+6:00
on 1996-10-26, moving from 00:29:59 back to 00:00:00.

```julia
julia> lzdt = LaxZonedDateTime(DateTime(1996, 10, 26, 0, 45), TimeZone("Asia/Colombo"))
1996-10-25T23:45:00+06:30

julia> floor(lzdt, Dates.Hour)
1996-10-26T00:00:00+06:00

julia> floor(lzdt, Dates.Day)
1996-10-26T00:00:00-AMB
```
"""
Base.floor(::LaxZonedDateTimes.LaxZonedDateTime, ::Union{Period, Type{Period}})

"""
    ceil(lzdt::ZonedDateTime, p::Period) -> LaxZonedDateTime
    ceil(lzdt::ZonedDateTime, p::Type{Period}) -> LaxZonedDateTime

Returns the nearest `LaxZonedDateTime` greater than or equal to `lzdt` at resolution `p`.
The result will be in the same time zone as `lzdt`.

For convenience, `p` may be a type instead of a value: `ceil(lzdt, Dates.Hour)` is a
shortcut for `ceil(lzdt, Dates.Hour(1))`.

`VariableTimeZone` transitions are handled as for `round`.

### Examples

The `America/Winnipeg` time zone transitioned from Central Standard Time (UTC-6:00) to
Central Daylight Time (UTC-5:00) on 2016-03-13, moving directly from 01:59:59 to 03:00:00.

```julia
julia> lzdt = LaxZonedDateTime(DateTime(2016, 3, 13, 1, 45), TimeZone("America/Winnipeg"))
2016-03-13T01:45:00-06:00

julia> ceil(lzdt, Dates.Day)
2016-03-14T00:00:00-05:00

julia> ceil(lzdt, Dates.Hour)
2016-03-13T03:00:00-05:00
```

The `Asia/Colombo` time zone revised the definition of Lanka Time from UTC+6:30 to UTC+6:00
on 1996-10-26, moving from 00:29:59 back to 00:00:00.

```julia
julia> lzdt = LaxZonedDateTime(DateTime(1996, 10, 25, 23, 45), TimeZone("Asia/Colombo"))
1996-10-25T23:45:00+06:30

julia> ceil(lzdt, Dates.Hour)
1996-10-26T00:00:00+06:30

julia> ceil(lzdt, Dates.Day)
1996-10-26T00:00:00-AMB
```
"""
Base.ceil(::LaxZonedDateTimes.LaxZonedDateTime, ::Union{Period, Type{Period}})

"""
    round(lzdt::ZonedDateTime, p::Period, [r::RoundingMode]) -> LaxZonedDateTime
    round(lzdt::ZonedDateTime, p::Type{Period}, [r::RoundingMode]) -> LaxZonedDateTime

Returns the `LaxZonedDateTime` nearest to `lzdt` at resolution `p`. The result will be in
the same time zone as `lzdt`. By default (`RoundNearestTiesUp`), ties (e.g., rounding 9:30
to the nearest hour) will be rounded up.

For convenience, `p` may be a type instead of a value: `round(lzdt, Dates.Hour)` is a
shortcut for `round(lzdt, Dates.Hour(1))`.

Valid rounding modes for `round(::TimeType, ::Period, ::RoundingMode)` are
`RoundNearestTiesUp` (default), `RoundDown` (`floor`), and `RoundUp` (`ceil`).

### `VariableTimeZone` Transitions

`LaxZonedDateTime`s, like `ZonedDateTime`s, are rounded in their local time zone (rather
than UTC). This ensures that rounding behaves as expected and is maximally meaningful.

If rounding were done in UTC, consider how rounding to the nearest day would be resolved for
non-UTC time zones: the result would be 00:00 UTC, which wouldn't be midnight local time.
Similarly, when rounding to the nearest hour in `Australia/Eucla (UTC+08:45)`, the result
wouldn't be on the hour in the local time zone.

When `p` is a `DatePeriod` rounding is done in the local time zone in a straightforward
fashion. When `p` is a `TimePeriod` the likelihood of encountering an ambiguous or
non-existent time (due to daylight saving time transitions) is increased. To resolve this
issue, rounding a `LaxZonedDateTime` with a `VariableTimeZone` to a `TimePeriod` uses the
`DateTime` value in the appropriate `FixedTimeZone`, then reconverts it to a
`LaxZonedDateTime` in the appropriate `VariableTimeZone` afterward.

While it might be tempting to suppose that both `DatePeriod`s and `TimePeriod`s should be
rounded in local time (allowing ambiguous or non-existent times to crop up where they may),
`LaxZonedDateTime` attempts to return the maintain consistency with `ZonedDateTime` wherever
possible. While this may be contentious, most parties can probably agree that 01:59 rounded
to the nearest hour on a "spring forward" day should return 03:00 (the next hour that
exists) instead of 02:00 (which is non-existent).

Rounding is not an entirely "safe" operation for `LaxZonedDateTime`s, as in some cases
historical transitions for some time zones (such as `Asia/Colombo`) occur at midnight. In
such cases rounding to a `DatePeriod` may still return a result marked as ambiguous or
non-existent. Unlike a `ZonedDateTime`, however, no exceptions will be thrown.

### Rounding Ambiguous or Non-Existent Times

Rounding a `LaxZonedDateTime` that is in an ambiguous or non-existent state to a
`DatePeriod` will function as expected (because most transitions do not occur on date
boundaries, this will typically return the `LaxZonedDatetime` to a valid state).

Rounding an ambiguous or non-existent `LaxZonedDateTime` to a `TimePeriod` will result in
a non-representable state that is unrecoverable.

When `round`ing a `LaxZonedDateTime` using `RoundNearestTiesUp`, `round` needs to decide
whether `floor` or `ceil` is "closer" to the original date. In cases where `floor` or `ceil` result
would return a result that is either ambiguous or non-existent, `round` is unable to
determine which result is to be preferred, and will return a non-representable
`LaxZonedDateTime`.

### Examples

The `America/Winnipeg` time zone transitioned from Central Standard Time (UTC-6:00) to
Central Daylight Time (UTC-5:00) on 2016-03-13, moving directly from 01:59:59 to 03:00:00.

```julia
julia> lzdt = LaxZonedDateTime(DateTime(2016, 3, 13, 1, 45), TimeZone("America/Winnipeg"))
2016-03-13T01:45:00-06:00

julia> round(lzdt, Dates.Hour)
2016-03-13T03:00:00-05:00

julia> round(lzdt, Dates.Day)
2016-03-13T00:00:00-06:00

julia> lzdt = LaxZonedDateTime(DateTime(2016, 3, 13, 2, 45), TimeZone("America/Winnipeg"))
2016-03-13T02:45:00-DNE

julia> round(lzdt, Dates.Hour)
INVALID
```

The `Asia/Colombo` time zone revised the definition of Lanka Time from UTC+6:30 to UTC+6:00
on 1996-10-26, moving from 00:29:59 back to 00:00:00.

```julia
julia> lzdt = LaxZonedDateTime(DateTime(1996, 10, 25, 23, 45), TimeZone("Asia/Colombo"))
1996-10-25T23:45:00+06:30

julia> round(lzdt, Dates.Hour)
1996-10-26T00:00:00+06:30

julia> round(lzdt, Dates.Day)
INVALID
```
"""     # Defined in base/dates/rounding.jl
Base.round(::LaxZonedDateTimes.LaxZonedDateTime, ::Union{Period, Type{Period}})