text:
# function drv(p,n=1,isfreq=false)
#   d=Categorical(p)
#   if n==1
#     return rand(d)
#   end
#   if isfreq
#     pn = round(p[1:end-1]*n)
#     pn = [pn,n-sum(pn)]
#
#     pnc = zeros(Int,length(p))
#     ps = zeros(Int,n)
#     for i in 1:n
#       while true
#         v = rand(d)
#         if pnc[v] < pn[v]
#           pnc[v] += 1
#           ps[i]=v
#           break
#         end
#       end
#     end
#     return ps
#   else
#     rand(d,n)
#   end
# end
#
# function randvec3(xr::Vector,yr::Vector,zr::Vector)
#   dx = Uniform(xr[1],xr[2])
#   dy = Uniform(yr[1],yr[2])
#   dz = Uniform(zr[1],zr[2])
#   Vec3(rand(dx),rand(dy),rand(dz))
# end
# randvec3(;xmin=0,xmax=1,ymin=0,ymax=1,zmin=0,zmax=1)=randvec3([xmin,xmax],[ymin,ymax],[zmin,zmax])
# function randvec3(absmin::Vector)
#   d = Uniform(-1,1)
#   av3 = zeros(3)
#   for i=1:3
#     while true
#       r= rand(d)
#       if abs(r) >= absmin[i]
#         av3[i]=r
#         break;
#       end
#     end
#   end
#   return convert(Vec3,av3)
# end
# randvec3(xabsmin,yabsmin,zabsmin)=randvec3([xabsmin,yabsmin,zabsmin])
# function randvec3(radmin;is2d=false)
#   innerabsmin = radmin*sin(pi/4)
#   absmin = fill(innerabsmin,3)
#   while true
#     cv3 = randvec3(absmin)
#     if is2d;cv3.z=0.0;end
#     if length(cv3) >= radmin
#       return cv3
#     end
#   end
# end

# function convert(::Type{DataFrame},ct::CoefTable)
#    df = convert(DataFrame,ct.mat)
#    names!(df,map(symbol,ct.colnms))
#    [DataFrame(coefname=ct.rownms) df]
#end

# function flcond(fl,fs...)
#   if length(fs)==0
#     fs=collect(keys(fl))
#   end
#   ex = :([{(fs[1],l1)} for l1=fl[fs[1]]][:])
#   for i=2:length(fs)
#     lx=symbol("l$i")
#     push!(ex.args[1].args,:($lx=fl[fs[$i]]))
#     push!(ex.args[1].args[1].args,:((fs[$i],$lx)))
#   end
#   eval(ex)
# end
function flcond(fl,f1)
    conds = [Any[(f1,l1)] for l1 in fl[f1]]
end
function flcond(fl,f1,f2)
    conds = [Any[(f1,l1),(f2,l2)] for l1 in fl[f1],l2 in fl[f2]][:]
end
function flcond(fl,f1,f2,f3)
    conds = [Any[(f1,l1),(f2,l2),(f3,l3)] for l1 in fl[f1],l2 in fl[f2],l3 in fl[f3]][:]
end
function subcond(conds,sc...)
    sci = map(i->issubset(sc,i),conds)
    return conds[sci]
end

"Find Condition with Factor=level"
function findcond(df::DataFrame;fl...)
    i = trues(size(df,1))
    for f in keys(fl)
        i .&= df[!,f].==fl[f]
    end
    return df[findall(i),:]
end
function findcond(df::DataFrame,cond::Vector{Any};roundingdigit=3)
    i = trues(size(df,1))
    condstr = ""
    for fl in cond
        f = fl[1]
        l = fl[2]
        i &= df[Symbol(f)].==l
        if typeof(l)<:Real
            lv=round(l,roundingdigit)
        else
            lv=l
        end
        condstr = "$condstr, $f=$lv"
    end
    return find(i),condstr[3:end]
end
function findcond(df::DataFrame,conds::Vector{Vector{Any}};roundingdigit=3)
    n = length(conds)
    is = Array(Vector{Int},n)
    ss = Array(String,n)
    for i in 1:n
        is[i],ss[i] = findcond(df,conds[i],roundingdigit=roundingdigit)
    end
    vi = map(i->!isempty(i),is)
    return is[vi],ss[vi],conds[vi]
end

"Check if `response` is significently different from `baseline` by `Wilcoxon Signed Rank Test`"
isresponsive(baseline,response;alpha=0.05) = pvalue(SignedRankTest(baseline,response)) < alpha
"Check if any `sub group of response` is significently different from `baseline` by `Wilcoxon Signed Rank Test`"
isresponsive(baseline,response,gi;alpha=0.05) = any(map(i->isresponsive(baseline[i],response[i],alpha=alpha),gi))
isresponsive(baseline::Vector,response::Matrix;alpha=0.05) = any(isresponsive.(baseline,response,alpha=alpha))
"Check if `std` of a spatial-temporal kernal within response time window significently different from that of the baseline time window"
function isresponsive(st::Matrix,bti::Vector;sdfactor=3)
    sd = dropdims(std(st,dims=1),dims=1)
    sdmaxt = argmax(sd); sdmax = sd[sdmaxt]
    bsdm = mean(sd[bti]);bsdsd=std(sd[bti])

    (;r=!(sdmaxt in bti) && sdmax > bsdm+sdfactor*bsdsd,sd=sdmax,d=sdmaxt)
end

"Check if any factors and their interactions significently modulate response using ANOVA"
function ismodulative(df;alpha=0.05,interact=true)
    xns = filter(i->i!=:Y,propertynames(df))
    foreach(i->df[!,i]=categorical(df[!,i]),xns)
    if interact
        f = term(:Y) ~ reduce(+,map(i->reduce(&,term.(i)),combinations(xns)))
    else
        f = term(:Y) ~ reduce(+,term.(xns))
    end
    lmr = fit(LinearModel,f,df,contrasts = Dict(x=>EffectsCoding() for x in xns))
    anovatype = length(xns) <= 1 ? 2 : 3
    any(Anova(lmr,anovatype = anovatype).p[1:end-1] .< alpha)
end

"""
Find levels for each factor and indices, repetition for each level
"""
flin(ctc::Dict)=flin(DataFrame(ctc))
function flin(ctc::DataFrame)
    isempty(intersect(propertynames(ctc), (:i,:n))) || error("i and n are reserved for condition test indices and repeats, shouldn't be used for factor name.")
    fl=OrderedDict{Symbol,DataFrame}()
    for f in propertynames(ctc)
        fl[f] = condin(ctc[:,[f]])
    end
    return fl
end

"""
In Condition Tests, find each unique condition, number of repeats and its indices
"""
condin(ctc::Dict)=condin(DataFrame(ctc))
function condin(ctc::DataFrame)
    isempty(intersect(propertynames(ctc), (:i,:n))) || error("i and n are reserved for condition test indices and repeats, shouldn't be used for factor name.")
    combine(groupby([ctc DataFrame(i=1:nrow(ctc))], propertynames(ctc), sort = true, skipmissing = true), :i => (x->[x]) => :i, nrow => :n)
end

"Get factors of conditions"
condfactor(cond)=setdiff(propertynames(cond),(:i,:n))

"Get `Final` factors of conditions"
function finalfactor(cond::DataFrame)
    fs = String.(condfactor(cond))
    filter!(f->endswith(f,"_Final") || ∉("$(f)_Final",fs),fs)
    Symbol.(fs)
end

"Print Condition in String"
function condstring(cond::DataFrameRow;factor=condfactor(cond))
    join(["$f=$(cond[f])" for f in factor],", ")
end
function condstring(cond::DataFrame;factor=condfactor(cond))
    [condstring(r,factor=factor) for r in eachrow(cond)]
end

"""
Group repeats of Conditions, get `Mean` and `SEM` of responses

1. rs: responses of each trial
2. ci: trial indices of repeats for each condition
"""
function condresponse(rs,ci)
    crs = [rs[i] for i in ci]
    DataFrame(m=mean.(crs),se=sem.(crs))
end
function condresponse(rs,cond::DataFrame;u=0,ug="SU")
    crs = [rs[i] for i in cond.i]
    [DataFrame(m=mean.(crs),se=sem.(crs),u=fill(u,length(crs)),ug=fill(ug,length(crs))) cond[:,condfactor(cond)]]
end
function condresponse(urs::Dict,cond::DataFrame)
    mapreduce(k->condresponse(urs[k],cond,u=k),append!,keys(urs))
end
function condresponse(urs::Dict,ctc::DataFrame,factors)
    vf = intersect(propertynames(ctc),factors)
    isempty(vf) && error("No Valid Factor Found.")
    condresponse(urs,condin(ctc[:,vf]))
end
"""
Mean image responses of each conditions

1. rs: responses of each epoch [height, width, nepoch]
2. ci: epoch indices of repeats for each condition
"""
function condimageresponse(rs,ci)
    cr = Array{Float64}(undef,size(rs)[1:2]...,length(ci))
    @inbounds for i in eachindex(ci)
        @views cr[:,:,i] = dropdims(mean(rs[:,:,ci[i]],dims=3),dims=3)
    end
    cr
end

"Condition Response in Factor Space"
function factorresponse(df;factors = setdiff(propertynames(df),[:m,:se,:u,:ug]),fl = flin(df[:,factors]),fa = OrderedDict(f=>fl[f][!,f] for f in keys(fl)))
    fm = missings(Float64, map(nrow,values(fl))...)
    fse = deepcopy(fm)
    for i in 1:nrow(df)
        idx = [findfirst(df[i:i,f].==fa[f]) for f in keys(fa)]
        fm[idx...] = df[i,:m]
        fse[idx...] = df[i,:se]
    end
    return fm,fse,fa
end
function factorresponse(unitspike,cond,condonoff,dataset,unitsync)
    fms=[];fses=[];factors = condfactor(cond)
    fl = flin(cond[:,factors]);fa = OrderedDict(f=>fl[f][!,f] for f in keys(fl))
    for u in eachindex(unitspike)
        rs = epochspiketrainresponse_ono(unitspike[u],ref2sync(condonoff,dataset,unitsync[u]),israte=true)
        df = condresponse(rs,cond)
        fm,fse,_ = factorresponse(df,factors=factors,fl=fl,fa=fa)
        push!(fms,fm);push!(fses,fse)
    end
    return fms,fses,fa
end

function setfln(fl::Dict,n::Int)
    fln=Dict()
    for f in keys(fl)
        for l in fl[f]
            fln[(f,l)] = n
        end
    end
    return fln
end
function setfln(fl::Dict,fn::Dict)
    fln=Dict()
    for f in keys(fl)
        for i=1:length(fl[f])
            fln[(f,fl[f][i])] = fn[f][i]
        end
    end
    return fln
end

function testfln(fln::Dict,minfln::Dict;showmsg::Bool=true)
    r=true
    for mkey in keys(minfln)
        if haskey(fln,mkey)
            if fln[mkey] < minfln[mkey]
                if showmsg;print("Level $(mkey[2]) of factor \"$(mkey[1])\" repeats $(fln[mkey]) < $(minfln[mkey]).\n");end
                r=false
            end
        else
            if showmsg;print("Level $(mkey[2]) of factor \"$(mkey[1])\" is missing.\n");end
            r=false
        end
    end
    return r
end
function testfln(ds::Vector,minfln::Dict;showmsg::Bool=true)
    vi = find(map(d->begin
    if showmsg;print("Testing factor/level of \"$(d["datafile"])\" ...\n");end
    testfln(d["fln"],minfln,showmsg=showmsg)
end,ds))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3247c742-65df-11ec-0eed-7d76c7786941
begin
	using Pkg; Pkg.activate(@__DIR__); Pkg.instantiate()
	Pkg.precompile()

	using PlutoUI
	using GLMNet
	using RDatasets
	using MLBase
	using Plots
	using DecisionTree
	using Distances
	using NearestNeighbors
	using Random
	using LinearAlgebra
	using DataStructures
	using LIBSVM
end;

# ╔═╡ ff42960a-f27f-41b8-8eee-14fb5b0454ba
PlutoUI.TableOfContents(aside=true, indent=true, depth=3)

# ╔═╡ b4c7b20b-ef19-4684-933b-7a00bfb69dce
html"""

"""

# ╔═╡ 73a1eed0-af5b-42b3-a4b8-a6f191f49daa
md"""
# 🔧 Classification

Put simply, classification is the task of predicting a label for a given observation. For example: you are given certain physical descriptions of an animal, and your taks is to classify them as either a dog or a cat. Here, we will classify iris flowers.

[YouTube Link](https://www.youtube.com/watch?v=OQRPeIQasdo)

As we will see later, we will use different classifiers and at the end of this notebook, we will compare them. We will define our accuracy function right now to get it out of the way. We will use a simple **accuracy** function that returns the ratio of the number of correctly classified observations (true positives) to the total number of predictions.
"""

# ╔═╡ 68c8634d-a052-4bd0-a101-5111dd9d4f68
findaccuracy(ŷ, y) = sum(ŷ .== y) / length(y)

# ╔═╡ 091524fe-ad86-45a3-ad25-370fe25a7a4e
md"""
## Load data

We will use Iris dataset from [RDatasets.jl](https://github.com/JuliaStats/RDatasets.jl) package. It is already stored as a `DataFrame`.
"""

# ╔═╡ 7a1a6481-eb71-4eb8-b47e-cf65c2c4cafa
iris = dataset("datasets", "iris")

# ╔═╡ 0bca5475-72f2-4189-b3d8-dd06cc951fc8
md"""
## Data wrangling

First, let's define our feature matrix `X`...
"""

# ╔═╡ e7f4e859-5e3c-439d-b069-46f17f2acf02
X = Matrix(iris[:, 1:4])

# ╔═╡ 0c1f6c24-ab74-494a-89d8-ba3fbda7fa24
md"""
And our label vector `labels`...
"""

# ╔═╡ b086ca45-6a60-4d4e-8387-de7e247812ed
labels = iris[:, 5]

# ╔═╡ c495a8f7-17d0-42f7-9adf-440d83d310c1
md"""
Since our dependent variable is categorical, let's encode it using the [MLBase.jl](https://github.com/JuliaStats/MLBase.jl) package...

**Note:** some algorithms require numeric labels rather than string labels.
"""

# ╔═╡ 35fad54f-e861-46e7-ab43-04a205c1991f
# mapping labels
irismap = labelmap(labels)

# ╔═╡ df121f29-a69d-42c9-b6a7-f129c4e94506
# enconding labels vector
y = labelencode(irismap, labels)

# ╔═╡ 1f39a7ed-7320-4f72-b4a0-fac44a71bb75
md"""
## Train/test split

In classification, we often want to use some of the data to fit a model (i.e. _training data_), and the rest of the data to validate (i.e. _testing data_). We will get this data ready now, so that we can easily use it in the rest of this notebook.

We will create the `train_test_split` function to perform train/test split. It will return the `y_train` indexes...

**Note:** the following `train_test_split` function performs a stratified random sampling by class.
"""

# ╔═╡ ecb46a59-1acf-48b9-b3c7-4c319388bc33
function train_test_split(y::Array{Int64,1}, train_ratio::Float64)
	classes = unique(y)
	y_train_idx = []
	for class in classes
		idxs = findall(y .== class)
		selected_idxs = randsubseq(idxs, train_ratio)
		push!(y_train_idx, selected_idxs...)
	end

	return y_train_idx
end

# ╔═╡ 8ca957e5-1bec-4498-aaca-86dccca1a49b
md"""
**Note:** `randsubseq(A, p)` function returns a vector consisting of a random subsequence of the given array `A`, where each element of `A` is included (in order) with independent probability `p`. Technically, this process is known as _Bernoulli sampling_ of `A`.
"""

# ╔═╡ 8868df3e-5933-4670-bbaf-eb8be100382b
# y train indexes
train_idxs = train_test_split(y, 0.7)

# ╔═╡ 27f97d52-c5c2-41a3-9faf-a8271a404b23
# y test indexes
test_idxs = setdiff(1:length(y), train_idxs)

# ╔═╡ 95a6bc1e-7986-4d86-9843-ceac5b668e4c
md"""
We will need one more function, and that is the function that will assign classes based on the predicted values when the predicted values are continuous.

**Example:** suppose that we have $\hat{y}_i = 0.95$. This is a continuous prediction and we may convert it into $1$ (i.e. _setosa_).
"""

# ╔═╡ 137c8d77-d6fd-44eb-937e-5978d201670d
assign_class(ŷᵢ) = argmin(abs.(ŷᵢ .- [1,2,3]))

# ╔═╡ 6317039e-1926-4143-a4d6-61f61b5e83e0
md"""
## Data modeling

We are going to build 7 different classification models:
1. Lasso
2. Ridge
3. Elastic Net
4. Decision Trees
5. Random Forests
6. Nearest Neighbour
7. Support Vector Machine

Then, we will compare model accuracies to figure out which models have the best performance.
"""

# ╔═╡ 1bd96a25-9542-432b-8e54-d6896a9cdade
md"""
### Lasso

**Least Absolute Shrinkage and Selection Operator (Lasso)** is a regression analysis method that performs both variable selection and regularization in order to enhance the prediction accuracy and interpretability of the resulting statistical model. It was originally introduced in Geophysics, and later by Robert Tibshirani who coined the term.

Lasso was originally formulated for linear regression models. This simple case reveals a substantial amount about the estimator. These include its relationship to Ridge and best subset selection and the connections between Lasso coefficient estimates and so-called _soft thresholding_. It also reveals that (like standard linear regression) the coefficient estimates do not need to be unique if covariates (independent variables) are collinear.

**Note:** check [here](https://en.wikipedia.org/wiki/Lasso_(statistics)) for more math details.

For Lasso, Ridge and Elastic Net methods, we are going to use the [GLMNet.jl](https://github.com/JuliaStats/GLMNet.jl) package.

We use the `glmnet` function, we will get several different solutions and $\lambda$ values (weighting parameter for the norm factor). Thus, we may use the `glmnetcv` function to get the _best_ $\lambda$ value (i.e. the $\lambda$ value with the lowest mean loss).

**Note:** to perform Lasso, we must set `alpha=1`, which is the default `alpha` value for the `glmnet` function.
"""

# ╔═╡ dd6889f6-18da-4218-bc48-b6589699498a
# Lasso solutions
lasso_sol = glmnet(X[train_idxs, :], y[train_idxs])

# ╔═╡ d03a03d1-fc97-4cb3-a852-5a3bc359a141
# performing cross validation
lasso_cv = glmnetcv(X[train_idxs, :], y[train_idxs])

# ╔═╡ d6008ed9-2a46-446a-b684-1d961c3ae20a
# best λ for Lasso
lasso_λ = lasso_sol.lambda[argmin(lasso_cv.meanloss)]

# ╔═╡ 039b7860-02fb-4f97-a1c2-9d20af989aac
# Lasso final model
lasso_f̂ = glmnet(X[train_idxs, :], y[train_idxs],lambda=[lasso_λ])

# ╔═╡ d0192106-f984-41e8-aa50-897481729dca
md"""
Now, we can predict labels for testing data...
"""

# ╔═╡ c106c3da-6444-4e64-9530-b6102766cd22
lasso_pred = GLMNet.predict(lasso_f̂, X[test_idxs,:])

# ╔═╡ 0da3fe69-f22d-4a29-9c88-c6b8e2e13d08
md"""
Note that predictions above are continuous rather than discrete. Thus, we must use the `assign_class` function to convert continuous predictions into discrete predictions...
"""

# ╔═╡ 1da8374e-aba5-4da3-bfe4-8f450408a20b
lasso_ŷ = assign_class.(lasso_pred)

# ╔═╡ 26ff22e5-e54e-4b85-b6d2-0db838cc1aca
md"""
Finally, we can compute Lasso accuracy using the `findaccuracy` function...
"""

# ╔═╡ 5d31a589-459e-4c11-8ea8-9840eb2cf8c0
lasso_acc = findaccuracy(lasso_ŷ, y[test_idxs])

# ╔═╡ 0040a4e5-13d3-4277-a865-e3f1bd328388
md"""
### Ridge

**Ridge** is a method of estimating the coefficients of multiple-regression models in scenarios where independent variables are highly correlated.

Ridge regression was developed as a possible solution to the imprecision of least square estimators when linear regression models have some multicollinear (highly correlated) independent variables—by creating a ridge regression estimator (RR). This provides a more precise ridge parameters estimate, as its variance and mean square estimator are often smaller than the least square estimators previously derived.

The main difference between Lasso and Ridge is that the former optimizes for L1-norm and Ridge optimizes for L2-norm.

**Note:** check [here](https://en.wikipedia.org/wiki/Ridge_regression) for more math details.

To perform Ridge, we are going to use the same functions. However, we must set `alpha = 0`.
"""

# ╔═╡ face4238-2511-4de5-8df3-2d25a668386e
# Ridge solutions
ridge_sol = glmnet(X[train_idxs, :], y[train_idxs], alpha=0)

# ╔═╡ faf82a4f-5acb-4d45-bf33-ded1dd82733c
# performing cross validation
ridge_cv = glmnetcv(X[train_idxs, :], y[train_idxs], alpha=0)

# ╔═╡ 103b74be-e907-41c5-ba08-3bd3ee4242e1
# best λ for Ridge
ridge_λ = ridge_sol.lambda[argmin(ridge_cv.meanloss)]

# ╔═╡ b16bfae9-5020-48fd-a796-f77d67a53f2a
# Ridge final model
ridge_f̂ = glmnet(X[train_idxs, :], y[train_idxs],lambda=[ridge_λ], alpha=0)

# ╔═╡ 4d175b94-43d1-497c-8022-42f544789e6d
md"""
Now, we can predict labels for testing data...
"""

# ╔═╡ e5e6130f-62de-4966-b164-9323b0b3bfb5
ridge_pred = GLMNet.predict(ridge_f̂, X[test_idxs,:])

# ╔═╡ a5a16631-0cd9-44cf-ae12-3f66b822f74b
md"""
Note that predictions above are also continuous rather than discrete. Thus, we must use the `assign_class` function to convert continuous predictions into discrete predictions...
"""

# ╔═╡ 887520af-9c11-4185-9f97-7da958332f76
ridge_ŷ = assign_class.(ridge_pred)

# ╔═╡ 6e37ae05-527f-4e5d-9207-383c0290f8d9
md"""
Finally, we can compute Ridge accuracy using the `findaccuracy` function...
"""

# ╔═╡ cf702060-5b49-451a-9075-3fb6f544e666
ridge_acc = findaccuracy(ridge_ŷ, y[test_idxs])

# ╔═╡ 94547be0-0467-4323-b5cd-f1efa9b98ac8
md"""
### Elastic Net

In the fitting of linear or logistic regression models, the **Elastic Net** is a regularized regression method that linearly combines the L1 and L2 penalties of the lasso and ridge methods.

**Note:** check [here](https://en.wikipedia.org/wiki/Elastic_net_regularization) for more math details.

To perform Elastic Net, we are going to use the same functions. However, we must set `alpha = 0.5`.
"""

# ╔═╡ ca590a9a-5909-4f1e-b02e-8af5da63439f
# Elastic Net solutions
enet_sol = glmnet(X[train_idxs, :], y[train_idxs], alpha=0.5)

# ╔═╡ d2e707e7-81ef-4b3d-8b8b-5ef65c8fc754
# performing cross validation
enet_cv = glmnetcv(X[train_idxs, :], y[train_idxs], alpha=0.5)

# ╔═╡ 9da0251a-ec4e-4b15-9527-9212af497e61
# best λ for Elastic Net
enet_λ = enet_sol.lambda[argmin(enet_cv.meanloss)]

# ╔═╡ eb72d5c2-15c4-41a9-81a2-c7c9f5f047ec
# Elastic Net final model
enet_f̂ = glmnet(X[train_idxs, :], y[train_idxs], alpha=0.5, lambda=[enet_λ])

# ╔═╡ 7a73633e-87d2-449b-8115-db6bf7e00266
md"""
Now, we can predict labels for testing data...
"""

# ╔═╡ 87906726-319e-454e-9299-77f6647081b8
enet_pred = GLMNet.predict(enet_f̂, X[test_idxs, :])

# ╔═╡ 5cdc7e8c-42e4-4f75-ab10-0a2b65d55809
md"""
Note that predictions above are also continuous rather than discrete. Thus, we must use the `assign_class` function to convert continuous predictions into discrete predictions...
"""

# ╔═╡ 14e764bb-94d9-4b1a-b2ce-1a39b930334a
enet_ŷ = assign_class.(enet_pred)

# ╔═╡ 46e9152e-7ee0-4420-8de6-331cef06b6dd
md"""
Finally, we can compute Elastic Net accuracy using the `findaccuracy` function...
"""

# ╔═╡ de00cd9a-009b-4991-aace-17ea06e41a1a
enet_acc = findaccuracy(enet_ŷ, y[test_idxs])

# ╔═╡ d999d681-d250-401e-87fe-963e625bfa56
md"""
### Decision Trees

**Decision Trees Learning** uses a decision tree (as a predictive model) to go from observations about an item (represented in the branches) to conclusions about the item's target value (represented in the leaves). Tree models where the target variable can take a discrete set of values are called classification trees; in these tree structures, leaves represent class labels and branches represent conjunctions of features that lead to those class labels.

**Note:** check [here](https://en.wikipedia.org/wiki/Decision_tree) for more math details.

We are going to use the [DecisionTree.jl](https://github.com/bensadeghi/DecisionTree.jl) package.
"""

# ╔═╡ 7936c964-03a4-4498-83ad-25a6cc7825a5
# defining Decision Tree model
dt_f̂ = DecisionTreeClassifier(max_depth=2)

# ╔═╡ 00258a2e-3c99-4810-879c-9b15489e52ae
# training algorithm
DecisionTree.fit!(dt_f̂, X[train_idxs,:], y[train_idxs])

# ╔═╡ 4c0e88ae-a7ca-45ca-9a23-fbd4cd752cbe
md"""
Now, we can predict labels for testing data...
"""

# ╔═╡ 02ef2267-6694-4006-b9c2-21241e7ff13a
dt_ŷ = DecisionTree.predict(dt_f̂, X[test_idxs, :])

# ╔═╡ da6959b5-29a0-42c6-8a55-3711b5f41074
md"""
Finally, we can compute Decision Tree accuracy using the `findaccuracy` function...
"""

# ╔═╡ bc55879d-d6b2-4562-a710-988968aa3ae4
dt_acc = findaccuracy(dt_ŷ, y[test_idxs])

# ╔═╡ 50a8aa34-fa23-4dbc-b69f-b897a54ad157
md"""
### Random Forests

**Random Forests** are an ensemble learning method for classification, regression and other tasks that operates by constructing a multitude of decision trees at training time. For classification tasks, the output of the random forest is the class selected by most trees. Random decision forests correct for decision trees habit of overfitting to their training set.

**Note:** check [here](https://en.wikipedia.org/wiki/Random_forest) for more math details.

The `RandomForestClassifier` function is available through the [DecisionTree.jl](https://github.com/bensadeghi/DecisionTree.jl) package as well.
"""

# ╔═╡ 71672bfa-558e-4ceb-88fb-aa61f615603e
# defining Random Forests model
rf_f̂ = RandomForestClassifier(n_trees=20)

# ╔═╡ 7d3e38ab-009a-445b-9e63-acc8e3ac4d8c
# training algorithm
DecisionTree.fit!(rf_f̂, X[train_idxs,:], y[train_idxs])

# ╔═╡ 18564590-adaa-4f6d-8180-0c4ece6a4398
md"""
Now, we can predict labels for testing data...
"""

# ╔═╡ b19ddb74-a299-48f7-9bec-20d078edafb4
rf_ŷ = DecisionTree.predict(rf_f̂, X[test_idxs,:])

# ╔═╡ a36aad87-8772-4ff6-9e63-34313346f37e
md"""
Finally, we can compute Decision Tree accuracy using the `findaccuracy` function...
"""

# ╔═╡ 04d476ca-07cb-4138-bd1b-679c223d7faf
rf_acc = findaccuracy(rf_ŷ, y[test_idxs])

# ╔═╡ 626d7979-0362-4a54-8f7d-1df787a68195
md"""
### Nearest Neighbour

The **Nearest Neighbour Search** algorithm aims to find the point in the tree that is nearest to a given input point. This search can be done efficiently by using the tree properties to quickly eliminate large portions of the search space.

**Note:** check [here](https://en.wikipedia.org/wiki/K-d_tree) for more math details.

We are going to use the [NearestNeighbors.jl](https://github.com/KristofferC/NearestNeighbors.jl) package. This is not straight forward, since we need to do some preprocessing first. Thus, let's build the KD tree.

**Note:** we need to transpose `X` matrix in order to use NearestNeighbors.jl package.
"""

# ╔═╡ f81dfbd1-b018-4098-a484-2619215abace
# building KD tree
kdtree = KDTree(X[train_idxs,:]')

# ╔═╡ 14c7f4ba-1e06-4fe5-9fb5-b0e9bbf4793d
md"""
Now, we will use the `knn` function, which performs a lookup of the `k` nearest neigbours to the `points` from the data in the `tree`. If `sortres = true`, the result is sorted such that the results are in the order of increasing distance to the point. It returns both `indices` and `distances` of the training data.
"""

# ╔═╡ 0dd179a1-d48d-42d2-85c6-62ba48af0ccb
# getting indices and distances
knn_idxs, knn_dists = knn(kdtree, X[test_idxs,:]', 5, true);

# ╔═╡ 588afb06-4f7b-458f-8037-6ef53b96b60b
# indices
knn_idxs

# ╔═╡ ff2fcea0-3588-4a79-9906-f3effa148e86
md"""
**Note:** this is a nested `Array`, where each element is an 1D `Array` with the **indexes** of the 5 nearest points. Each element represents a test observation.
"""

# ╔═╡ 75d5de77-08db-4bdd-b418-fe9fa9e6353d
# distances
knn_dists

# ╔═╡ 4644848e-c9ae-4742-9ce3-4747a85b1fd2
md"""
**Note:** this is also a nested `Array`, where each element is an 1D `Array` with the **distances** of the 5 nearest points. Each element represents a test observation. Each element is in ascending distance order.
"""

# ╔═╡ fa554c22-2156-43a6-8014-cbd98914f27d
# 5 nearest train observation labels to each test observation
nn_y_train = y[train_idxs][hcat(knn_idxs...)]

# ╔═╡ be1405be-750e-4a74-895a-ed56beb0cd36
# counting the number of labels per test observation
nn_possible_labels = map(i -> counter(nn_y_train[:, i]), 1:size(nn_y_train, 2))

# ╔═╡ 45625714-3ff0-4b46-8087-14ecbe672d08
md"""
Now, we can predict labels for testing data. Thus, we need to select the most frequent label per each test observation.
"""

# ╔═╡ ed93cb77-9eee-491a-835c-2952bd612946
# selecting the most frequent label per each test observation
nn_ŷ = map(i -> parse(Int,
					  string(string(argmax(nn_possible_labels[i])))),
		   1:size(nn_y_train, 2))

# ╔═╡ ec8854ea-542a-4e15-8533-ec6225cafa72
md"""
Finally, we can compute Nearest Neighbour accuracy using the `findaccuracy` function...
"""

# ╔═╡ 1ee8b7ec-3542-40ae-945c-3bf489266985
nn_acc = findaccuracy(nn_ŷ, y[test_idxs])

# ╔═╡ 112a337f-11ab-4074-88bd-17a887562d43
md"""
### Support Vector Machine

**Support Vector Machines (SVM)** are supervised learning models with associated learning algorithms that analyze data for classification and regression analysis.

Given a set of training examples, each marked as belonging to one of two categories, an SVM training algorithm builds a model that assigns new examples to one category or the other, making it a non-probabilistic binary linear classifier. SVM maps training examples to points in space so as to maximise the width of the gap between the two categories. New examples are then mapped into that same space and predicted to belong to a category based on which side of the gap they fall.

**Note:** check [here](https://en.wikipedia.org/wiki/Support-vector_machine) for more math details.

Here, we will use the [LIBSVM.jl](https://github.com/JuliaML/LIBSVM.jl) package. We need to transpose `X` matrix in order to use this package.
"""

# ╔═╡ 428463aa-0cc7-4008-95c7-8d1b5a28f2de
# training SVM model
svm_f̂ = svmtrain(X[train_idxs,:]', y[train_idxs])

# ╔═╡ fb93f2ce-36d2-48f0-bafd-5ca18afadfd7
md"""
Now, we can predict labels for testing data. At the same time, we will also get an `Array` with decision values...
"""

# ╔═╡ f5bb8271-4239-4c4e-956d-9714a38d2c25
svm_ŷ, decision_values = svmpredict(svm_f̂, X[test_idxs,:]');

# ╔═╡ 4a3d2073-f465-4fe2-9c4d-4ceb810a2285
# SVM predictions
svm_ŷ

# ╔═╡ ead6066e-6257-4659-84a6-98ea79993ea7
md"""
Finally, we can compute Support Vector Machine accuracy using the `findaccuracy` function...
"""

# ╔═╡ 92c46bcc-b3fa-4704-85cf-74782b9e6bef
svm_acc = findaccuracy(svm_ŷ, y[test_idxs])

# ╔═╡ 0f44eabe-c806-4154-b54e-65f8c678ef47
md"""
## Comparing models

Now, we can compare the accuracies of each model...
"""

# ╔═╡ e1fdde22-08ea-4ae4-82a6-b7adbb6e250a
begin
	models = ["Lasso", "Ridge", "Elastic Net", "DT", "RF", "NN", "SVM"]

	accs = [lasso_acc, ridge_acc, enet_acc, dt_acc, rf_acc, nn_acc, svm_acc]

	hcat(models, accs)
end

# ╔═╡ Cell order:
# ╟─3247c742-65df-11ec-0eed-7d76c7786941
# ╟─ff42960a-f27f-41b8-8eee-14fb5b0454ba
# ╟─b4c7b20b-ef19-4684-933b-7a00bfb69dce
# ╟─73a1eed0-af5b-42b3-a4b8-a6f191f49daa
# ╠═68c8634d-a052-4bd0-a101-5111dd9d4f68
# ╟─091524fe-ad86-45a3-ad25-370fe25a7a4e
# ╠═7a1a6481-eb71-4eb8-b47e-cf65c2c4cafa
# ╟─0bca5475-72f2-4189-b3d8-dd06cc951fc8
# ╠═e7f4e859-5e3c-439d-b069-46f17f2acf02
# ╟─0c1f6c24-ab74-494a-89d8-ba3fbda7fa24
# ╠═b086ca45-6a60-4d4e-8387-de7e247812ed
# ╟─c495a8f7-17d0-42f7-9adf-440d83d310c1
# ╠═35fad54f-e861-46e7-ab43-04a205c1991f
# ╠═df121f29-a69d-42c9-b6a7-f129c4e94506
# ╟─1f39a7ed-7320-4f72-b4a0-fac44a71bb75
# ╠═ecb46a59-1acf-48b9-b3c7-4c319388bc33
# ╟─8ca957e5-1bec-4498-aaca-86dccca1a49b
# ╠═8868df3e-5933-4670-bbaf-eb8be100382b
# ╠═27f97d52-c5c2-41a3-9faf-a8271a404b23
# ╟─95a6bc1e-7986-4d86-9843-ceac5b668e4c
# ╠═137c8d77-d6fd-44eb-937e-5978d201670d
# ╟─6317039e-1926-4143-a4d6-61f61b5e83e0
# ╟─1bd96a25-9542-432b-8e54-d6896a9cdade
# ╠═dd6889f6-18da-4218-bc48-b6589699498a
# ╠═d03a03d1-fc97-4cb3-a852-5a3bc359a141
# ╠═d6008ed9-2a46-446a-b684-1d961c3ae20a
# ╠═039b7860-02fb-4f97-a1c2-9d20af989aac
# ╟─d0192106-f984-41e8-aa50-897481729dca
# ╠═c106c3da-6444-4e64-9530-b6102766cd22
# ╟─0da3fe69-f22d-4a29-9c88-c6b8e2e13d08
# ╠═1da8374e-aba5-4da3-bfe4-8f450408a20b
# ╟─26ff22e5-e54e-4b85-b6d2-0db838cc1aca
# ╠═5d31a589-459e-4c11-8ea8-9840eb2cf8c0
# ╟─0040a4e5-13d3-4277-a865-e3f1bd328388
# ╠═face4238-2511-4de5-8df3-2d25a668386e
# ╠═faf82a4f-5acb-4d45-bf33-ded1dd82733c
# ╠═103b74be-e907-41c5-ba08-3bd3ee4242e1
# ╠═b16bfae9-5020-48fd-a796-f77d67a53f2a
# ╟─4d175b94-43d1-497c-8022-42f544789e6d
# ╠═e5e6130f-62de-4966-b164-9323b0b3bfb5
# ╟─a5a16631-0cd9-44cf-ae12-3f66b822f74b
# ╠═887520af-9c11-4185-9f97-7da958332f76
# ╟─6e37ae05-527f-4e5d-9207-383c0290f8d9
# ╠═cf702060-5b49-451a-9075-3fb6f544e666
# ╟─94547be0-0467-4323-b5cd-f1efa9b98ac8
# ╠═ca590a9a-5909-4f1e-b02e-8af5da63439f
# ╠═d2e707e7-81ef-4b3d-8b8b-5ef65c8fc754
# ╠═9da0251a-ec4e-4b15-9527-9212af497e61
# ╠═eb72d5c2-15c4-41a9-81a2-c7c9f5f047ec
# ╟─7a73633e-87d2-449b-8115-db6bf7e00266
# ╠═87906726-319e-454e-9299-77f6647081b8
# ╟─5cdc7e8c-42e4-4f75-ab10-0a2b65d55809
# ╠═14e764bb-94d9-4b1a-b2ce-1a39b930334a
# ╟─46e9152e-7ee0-4420-8de6-331cef06b6dd
# ╠═de00cd9a-009b-4991-aace-17ea06e41a1a
# ╟─d999d681-d250-401e-87fe-963e625bfa56
# ╠═7936c964-03a4-4498-83ad-25a6cc7825a5
# ╠═00258a2e-3c99-4810-879c-9b15489e52ae
# ╟─4c0e88ae-a7ca-45ca-9a23-fbd4cd752cbe
# ╠═02ef2267-6694-4006-b9c2-21241e7ff13a
# ╟─da6959b5-29a0-42c6-8a55-3711b5f41074
# ╠═bc55879d-d6b2-4562-a710-988968aa3ae4
# ╟─50a8aa34-fa23-4dbc-b69f-b897a54ad157
# ╠═71672bfa-558e-4ceb-88fb-aa61f615603e
# ╠═7d3e38ab-009a-445b-9e63-acc8e3ac4d8c
# ╟─18564590-adaa-4f6d-8180-0c4ece6a4398
# ╠═b19ddb74-a299-48f7-9bec-20d078edafb4
# ╟─a36aad87-8772-4ff6-9e63-34313346f37e
# ╠═04d476ca-07cb-4138-bd1b-679c223d7faf
# ╟─626d7979-0362-4a54-8f7d-1df787a68195
# ╠═f81dfbd1-b018-4098-a484-2619215abace
# ╟─14c7f4ba-1e06-4fe5-9fb5-b0e9bbf4793d
# ╠═0dd179a1-d48d-42d2-85c6-62ba48af0ccb
# ╠═588afb06-4f7b-458f-8037-6ef53b96b60b
# ╟─ff2fcea0-3588-4a79-9906-f3effa148e86
# ╠═75d5de77-08db-4bdd-b418-fe9fa9e6353d
# ╟─4644848e-c9ae-4742-9ce3-4747a85b1fd2
# ╠═fa554c22-2156-43a6-8014-cbd98914f27d
# ╠═be1405be-750e-4a74-895a-ed56beb0cd36
# ╟─45625714-3ff0-4b46-8087-14ecbe672d08
# ╠═ed93cb77-9eee-491a-835c-2952bd612946
# ╟─ec8854ea-542a-4e15-8533-ec6225cafa72
# ╠═1ee8b7ec-3542-40ae-945c-3bf489266985
# ╟─112a337f-11ab-4074-88bd-17a887562d43
# ╠═428463aa-0cc7-4008-95c7-8d1b5a28f2de
# ╟─fb93f2ce-36d2-48f0-bafd-5ca18afadfd7
# ╠═f5bb8271-4239-4c4e-956d-9714a38d2c25
# ╠═4a3d2073-f465-4fe2-9c4d-4ceb810a2285
# ╟─ead6066e-6257-4659-84a6-98ea79993ea7
# ╠═92c46bcc-b3fa-4704-85cf-74782b9e6bef
# ╟─0f44eabe-c806-4154-b54e-65f8c678ef47
# ╠═e1fdde22-08ea-4ae4-82a6-b7adbb6e250a

ext: jl

max_stars_repo_name: fnaghetini/DataScience

max_issues_repo_name: fnaghetini/DataScience

max_forks_repo_name: fnaghetini/DataScience


text:
            # P is the precision used for this value
type ArfFloat{P}  <: Real
  mid_exp::Int # fmpz
  mid_size::UInt # mp_size_t
  mid_d1::Int # mantissa_struct
  mid_d2::Int
end

precision{P}(x::ArfFloat{P}) = P
precision{P}(::Type{ArfFloat{P}}) = P

precision(::Type{ArfFloat}) = ArbFloatPrecision[1]
setprecision(::Type{ArfFloat}, x::Int) = setprecision(ArbFloat, x)

# a type specific hash function helps the type to 'just work'
const hash_arffloat_lo = (UInt === UInt64) ? 0x37e642589da3416a : 0x5d46a6b4
const hash_0_arffloat_lo = hash(zero(UInt), hash_arffloat_lo)
hash{P}(z::ArfFloat{P}, h::UInt) =
    hash(reinterpret(UInt,z.mid_d1)$z.mid_exp,
         (h $ hash(reinterpret(UInt,z.mid_d2)$(~reinterpret(UInt,P)), hash_arffloat_lo) $ hash_0_arffloat_lo))

function clearArfFloat{P}(x::ArfFloat{P})
     ccall(@libarb(arf_clear), Void, (Ptr{ArfFloat{P}},), &x)
end

function initializer{P}(::Type{ArfFloat{P}})
    z = ArfFloat{P}(0,0,0,0)
    ccall(@libarb(arf_init), Void, (Ptr{ArfFloat{P}},), &z)
    finalizer(z, clearArfFloat)
    z
end

midpoint{P}(x::ArfFloat{P}) = x

radius{P}(x::ArfFloat{P}) = zero(ArfFloat{P})

function frexp{P}(x::ArfFloat{P})
   mantissa = initializer(ArfFloat{P})
   exponent = zero(Int64)
   ccall(@libarb(arf_frexp), Void, (Ptr{ArfFloat{P}}, Int64, Ptr{ArfFloat{P}}), &mantissa, exponent, &x)
   mantissa, exponent
end


text:
using RanDAGen

dag = DAG(0:20, 2:4, 0:2)

# depth:  20
# width_min:  2       # [Sept 5]: 1, 2, 3
# width_max:  4       # [Sept 5]: 3, 4, 5
# ccr:  20            # [Sept 5]: 2, 20, 200
# comp_mu:    2       # computation workload
# comp_sigma: 0.5
# link_comm_sigma:  0.5   # [Sept 5] 0.1,0.5,0.9  #link_comm is in normal distribution. It is the data transmitted on links.
# deg_mu: 4           # [Sept 5]: 2, 4, 6
# deg_sigma: 0.5
# total_num: 5
# struct _Task {
#     int tag;
#     double cost;
#     int data_size;
#     double alpha;
#     int nb_children;
#     Task *children;
#     double *comm_costs;
#     int *transfer_tags;
#     complexity_t complexity;
#   };

#   struct _DAG {
#     int nb_levels;
#     int *nb_tasks_per_level;
#     Task **levels;
#   };


text:
const GroupIdx = UInt8

struct AgeCouplingParams
  coupling::CouplingSampler{Float64, PersonIdx, GroupIdx}
  source_weighting::Vector{Float64}
end

agegroup(age_thresholds::AbstractArray{T} where T<:Real, age::Real) = GroupIdx(searchsortedlast(age_thresholds, age))

function AgeCouplingParams(
  ages::AbstractVector{T} where T<:Real,
  genders::Union{Nothing, AbstractVector{Bool}},
  age_thresholds::AbstractVector{T} where T<:Real,
  coupling_weights::Matrix{Float64},
  contact_normalization::Union{Nothing, Real})

  @assert age_thresholds[1] == 0
  num_groups = length(age_thresholds)
  @assert (num_groups,num_groups) == size(coupling_weights)

  group_ids = agegroup.((age_thresholds,), ages)
  @assert minimum(group_ids) > 0
  if genders !== nothing
    @assert length(ages) == length(genders)
    group_ids .= group_ids .*2 .+ genders .- 1
  end
  @assert maximum(group_ids) <= num_groups <= typemax(GroupIdx) "minmax: $(extrema(group_ids)), num=$num_groups, typemax=$(typemax(GroupIdx))"

  source_weighting = Vector{Float64}(undef, size(coupling_weights, 2))
  sum!(source_weighting', coupling_weights)

  coupling_sampler = CouplingSampler(group_ids, coupling_weights, PersonIdx)

  if nothing !== contact_normalization # normalize such that expected number of contacts is equal to the given value
    group_sizes = groupsizes(coupling_sampler.grouping)
    mean_contacts = source_weighting' * group_sizes / sum(group_sizes)
    source_weighting .*= contact_normalization / mean_contacts
  end

  AgeCouplingParams(
    coupling_sampler,
    source_weighting
  )
end

sourceweightof(p::AgeCouplingParams, source_id::Integer) = p.source_weighting[groupof(p.coupling, source_id)]


text:
### SMMALA holds the fields of the SMMALA sampler
### These fields represent the initial user-defined state of the sampler
### These sampler fields are copied to the corresponding heap type fields, where the latter can be tuned

immutable SMMALA <: LMCSampler
  driftstep::Float64

  function SMMALA(ds::Float64)
    @assert ds > 0 "Drift step is not positive."
    new(ds)
  end
end

SMMALA(; driftstep::Float64=1.) = SMMALA(driftstep)

### SMMALAHeap type holds the internal state ("local variables") of the SMMALA sampler

type SMMALAHeap <: MCHeap{MCTensorSample}
  instate::MCState{MCTensorSample} # Monte Carlo state used internally by the sampler
  outstate::MCState{MCTensorSample} # Monte Carlo state outputted by the sampler
  tune::MCTune
  count::Int # Current number of iterations
  driftstep::Float64 # Drift stepsize for the current Monte Carlo iteration (possibly tuned)
  smean::Vector{Float64}
  ratio::Float64
  pnewgivenold::Float64
  poldgivennew::Float64
  current_invtensor::Matrix{Float64}
  successive_invtensor::Matrix{Float64}
  current_termone::Vector{Float64}
  successive_termone::Vector{Float64}
  cholinvtensor::Matrix{Float64}
end

SMMALAHeap() =
  SMMALAHeap(MCState(MCTensorSample(), MCTensorSample()), MCState(MCTensorSample(), MCTensorSample()), VanillaMCTune(),
  0, NaN, Float64[], NaN, NaN, NaN, Array(Float64, 0, 0), Array(Float64, 0, 0), Float64[], Float64[],
  Array(Float64, 0, 0))

SMMALAHeap(l::Int, t::MCTune=VanillaMCTune()) =
  SMMALAHeap(MCState(MCTensorSample(l), MCTensorSample(l)), MCState(MCTensorSample(l), MCTensorSample(l)), t, 0, NaN,
    fill(NaN, l), NaN, NaN, NaN, fill(NaN, l, l), fill(NaN, l, l), fill(NaN, l), fill(NaN, l),
    Array(Float64, 0, 0))

### Initialize SMMALA sampler

function initialize_heap(m::MCModel, s::SMMALA, r::MCRunner, t::MCTuner)
  @assert hasgradient(m) "SMMALA sampler requires model with gradient function."
  @assert hastensor(m) "SMMALA sampler requires model with tensor function."
  heap::SMMALAHeap = SMMALAHeap(m.size)

  heap.instate.current = MCTensorSample(copy(m.init))
  tensorlogtargetall!(heap.instate.current, m.evalallt)
  @assert isfinite(heap.instate.current.logtarget) "Initial values out of model support."

  heap.current_invtensor = inv(heap.instate.current.tensorlogtarget)
  heap.current_termone = heap.current_invtensor*heap.instate.current.gradlogtarget

  if isa(t, VanillaMCTuner)
    heap.tune = VanillaMCTune()
  elseif isa(t, EmpiricalMCTuner)
    heap.tune = EmpiricalMCTune(s.driftstep)
  end

  heap.count = 1

  heap
end

function reset!(heap::SMMALAHeap, x::Vector{Float64}, m::MCModel)
  heap.instate.current = MCTensorSample(copy(x))
  tensorlogtargetall!(heap.instate.current, m.evalallt)
end

function initialize_task!(heap::SMMALAHeap, m::MCModel, s::SMMALA, r::MCRunner, t::MCTuner)
  # Hook inside Task to allow remote resetting
  task_local_storage(:reset, (x::Vector{Float64})->reset!(heap, x, m))

  while true
    iterate!(heap, m, s, r, t, produce)
  end
end

### Perform iteration for SMMALA sampler

function iterate!(heap::SMMALAHeap, m::MCModel, s::SMMALA, r::MCRunner, t::MCTuner, send::Function)
  if isa(t, VanillaMCTuner)
    heap.driftstep = s.driftstep
    if t.verbose
      heap.tune.proposed += 1
    end
  elseif isa(t, EmpiricalMCTuner)
    heap.driftstep = heap.tune.step
    heap.tune.proposed += 1
  end

  # Calculate the drift term
  heap.smean = heap.instate.current.sample+(heap.driftstep/2)*heap.current_termone

  # Calculate proposed parameters
  heap.cholinvtensor = chol(heap.driftstep*heap.current_invtensor)
  heap.instate.successive = MCTensorSample(heap.smean+heap.cholinvtensor'*randn(m.size))

  # Update model based on the proposed parameters
  tensorlogtargetall!(heap.instate.successive, m.evalallt)

  heap.pnewgivenold = (-sum(log(diag(heap.cholinvtensor)))
    -(0.5*(heap.smean-heap.instate.successive.sample)'*(heap.instate.current.tensorlogtarget/heap.driftstep)
    *(heap.smean-heap.instate.successive.sample))[1])

  heap.successive_invtensor = inv(heap.instate.successive.tensorlogtarget)
  heap.successive_termone = heap.successive_invtensor*heap.instate.successive.gradlogtarget

  # Calculate the drift term
  heap.smean = heap.instate.successive.sample+(heap.driftstep/2)*heap.successive_termone

  heap.poldgivennew = (-sum(log(diag(chol(heap.driftstep*eye(m.size)*heap.successive_invtensor))))
    -(0.5*(heap.smean-heap.instate.current.sample)'*(heap.instate.successive.tensorlogtarget/heap.driftstep)
    *(heap.smean-heap.instate.current.sample))[1])

  heap.ratio = heap.instate.successive.logtarget+heap.poldgivennew-heap.instate.current.logtarget-heap.pnewgivenold

  if heap.ratio > 0 || (heap.ratio > log(rand()))
    heap.outstate = MCState(heap.instate.successive, heap.instate.current, Dict{Any, Any}("accept" => true))
    heap.instate.current = deepcopy(heap.instate.successive)
    heap.current_invtensor, heap.current_termone = copy(heap.successive_invtensor), copy(heap.successive_termone)

    if isa(t, VanillaMCTuner) && t.verbose
      heap.tune.accepted += 1
    elseif isa(t, EmpiricalMCTuner)
      heap.tune.accepted += 1
    end
  else
    heap.outstate = MCState(heap.instate.current, heap.instate.current, Dict{Any, Any}("accept" => false))
  end

  if isa(t, VanillaMCTuner) && t.verbose && heap.count <= r.burnin && mod(heap.count, t.period) == 0
    rate!(heap.tune)
    println("Burnin iteration $(heap.count) of $(r.burnin): ", round(100*heap.tune.rate, 2), " % acceptance rate")
  elseif isa(t, EmpiricalMCTuner) && heap.count <= r.burnin && mod(heap.count, t.period) == 0
    adapt!(heap.tune, t)
    reset!(heap.tune)
    if t.verbose
      println("Burnin iteration $(heap.count) of $(r.burnin): ", round(100*heap.tune.rate, 2), " % acceptance rate")
    end
  end

  heap.count += 1

  send(heap.outstate)
end


text:
module Competition

using RecordedArrays
using Random
using Test
using EvolutionaryModelingTools
using EvolutionaryModelingTools: sample

# parameters
clock = ContinuousClock(100.0)
r = 0.5   # growth rate
K = 100   # carrying capacity
μ = 0.01  # mutation rate

# state
v = recorded(DynamicEntry, clock, [100])
m = recorded(StaticEntry, clock, ones(1, 1))

# simulate with this package
@reaction growth begin
    r * (1 - μ) * v
    v[ind] += 1
end

@reaction mutation begin
    @. r * μ * v
    begin
        push!(v, 1)
        n = length(v)
        resize!(m, (n, n))
        # rand competition matrix
        m[:, n] = rand(rng, n)
        m[n, 1:n-1] = rand(rng, 1, n - 1)
    end
end

@reaction competition begin
    begin
        c = r / K # baseline competition coefficient
        @. c * v * m * v'
    end
    begin
        i = ind[1]
        v[i] -= 1
        if v[i] == 0 # check extinction
            deleteat!(v, i)
            resize!(m, (not(i), not(i)))
        end
    end
end

ps = (; r, K, μ, v, m)
rs = (growth, mutation, competition)

Random.seed!(1)
ps_updated = gillespie(clock, ps, rs)[1]

# simulate manually
Random.seed!(1)
init!(clock)

for _ in clock
    growth_rate = r * (1 - μ) * v
    mutation_rate = r * μ * v
    c = r / K # baseline competition coefficient
    competition_rate = @. c * v * m * v'
    rates_sum = (sum(growth_rate), sum(mutation_rate), sum(competition_rate))
    rates_sum_acc = accumulate(+, rates_sum)
    rates_sum_sum = rates_sum_acc[end]
    τ = -log(rand()) / rates_sum_sum
    increase!(clock, τ)
    rn = rand() * rates_sum_sum
    ind = searchsortedfirst(collect(rates_sum_acc), rn)
    if ind == 1
        i = sample(growth_rate, rn)
        v[i] += 1
    elseif ind == 2
        i = sample(mutation_rate, rn - rates_sum_acc[1])
        push!(v, 1)
        n = length(v)
        resize!(m, (n, n))
        # rand competition matrix
        m[:, n] = rand(n)
        m[n, 1:n-1] = rand(1, n - 1)
    elseif ind == 3
        i = sample(competition_rate, rn - rates_sum_acc[2])[1]
        v[i] -= 1
        if v[i] == 0 # check extinction
            deleteat!(v, i)
            resize!(m, (not(i), not(i)))
        end
    end
    length(v) == 0 && break
end

@test ps_updated.v == v
@test ps_updated.m == m

@test all(map(==, getentries(ps_updated.v), getentries(v)))
@test all(map(==, values(getentries(ps_updated.m)), values(getentries(m))))

end


text:
module IOExtension

export openmultiple

"""
    openmultiple(f::Function, files::Tuple{Vararg}, args...; kwargs....)

Apply the function `f` to various files, opened by broadcasting `open`
on a the `files` tuple and the additional arguments.
The resulting file descriptors are closed upon completion.
"""
function openmultiple(f::Function, files::Tuple{Vararg}, args...; kwargs...)
    io = open.(files, args...; kwargs...)
    try
        f(io...)
    finally
        close.(io)
    end
end

end # module


text:
using TestFunctionRunner

let env_test_cuda = lowercase(get(ENV, "UNIONARRAYS_JL_TEST_CUDA", "auto")),
    test_cuda = if env_test_cuda == "auto"
        try
            import CUDA
            true
        catch
            false
        end
    else
        env_test_cuda == "true"
    end

    if test_cuda
        TestFunctionRunner.@run(packages = ["UnionArraysCUDATests"])

        using FoldsCUDATests
        FoldsCUDATests.runtests_unionarrays()
    else
        TestFunctionRunner.@run(packages = ["UnionArraysTests"])
    end
end


text:
export flattened_tup_chain

flattened_tup_chain(::Type{NamedTuple{(), Tuple{}}}; prefix = (Symbol(),)) = ()
flattened_tup_chain(::Type{T}; prefix = (Symbol(),)) where {T <: Real} = (prefix,)
flattened_tup_chain(::Type{T}; prefix = (Symbol(),)) where {T <: SVector} = (prefix,)
flattened_tup_chain(::Type{T}; prefix = (Symbol(),)) where {T <: SHermitianCompact} = (prefix,)

"""
    flattened_tup_chain(::Type{T}) where {T <: Union{NamedTuple,NTuple}}

"""
function flattened_tup_chain(::Type{T}; prefix = (Symbol(),)) where {T <: Union{NamedTuple,NTuple}}
    map(1:fieldcount(T)) do i
        Ti = fieldtype(T, i)
        name = fieldname(T, i)
        sname = name isa Int ? name : Symbol(name)
        flattened_tup_chain(
            Ti,
            prefix = prefix == (Symbol(),) ? (sname,) : (prefix..., sname),
        )
    end |>
    Iterators.flatten |>
    collect
end


text:
"""
    ∇xJ(f, args...)

Returns gradients of the loss function (f) w.r.t. your parameters (args).
"""
function ∇xJ(f, args...)
    args = collect(args)
    args = map(a->isa(a, FloatArray) ? param(a) : a, args)
    J = @diff f(args...)
    ∇ = map(arg->grad(J, arg), args)
    ∇ = filter(i->isa(i, FloatArray), ∇)
end

"""
This function implements FGSM (Fast Gradient Sign Method)
 (arxiv.org/abs/1412.6572). FGSM is a white-box attack which means the attacker
 has full access to the model  including its parameters. In this attack, we
 simply take advantage of the  gradients of the loss with respect to input data.
 We first apply sign function to the input gradients, scale it, and then add
 those gradients to the input  data to get the perturbed data which will
 maximize the loss value. See lenet-fgsm notebook for example usage.

    FGSM(x::FloatArray, ∇x::FloatArray, ϵ;
         minval=0.0, maxval=1.0, targeted=false)

Applies FGSM for given input array (x), the gradient of your model's loss w.r.t.
input array (∇x) and ϵ. minval and maxval stands for clipping. If you want to
perform a targeted attack, use targeted=true.

    FGSM(f, ϵ, args...; minval=0.0, maxval=1.0, targeted=false)

Applies FGSM for given loss function (f), ϵ and arguments for the loss function.
"""
function FGSM(x::FloatArray, ∇x::FloatArray, ϵ;
              minval=0.0, maxval=1.0, targeted=false)
    T = eltype(x)
    sgn∇x = sign.(∇x)
    x̂ = x + T(ϵ) * sgn∇x * (-1.)^targeted
    x̂ = min.(T(maxval), x̂)
    x̂ = max.(T(minval), x̂)
    return x̂
end

function FGSM(f, ϵ, args...; minval=0.0, maxval=1.0, targeted=false)
    args = collect(args)
    ∇x  = ∇xJ(f, args...)
    x = filter(i->isa(i, FloatArray), args)
    x̂ = [FGSM(xi, ∇xi, ϵ; minval=minval, maxval=maxval, targeted=targeted)
         for (xi, ∇xi) in zip(x, ∇x)]
end

"""
    IterativeFGSM(f, ϵ, T::Int, args...;
                  minval=0.0, maxval=1.0, targeted=false)

Applies Iterative FGSM attack. Simply collects and returns adversarial
attack for all timesteps. f is a closure function which calls your models
loss function and only accept parameters you want to get gradients.
"""
function IterativeFGSM(f, ϵ, T::Int, args...; minval=0.0, maxval=1.0)
    x̂ = collect(args)
    history = []
    for t = 1:T
        x̂ = FGSM(f, ϵ/T, x̂...; minval=minval, maxval=maxval)
        push!(history, x̂)
    end
    return history
end

function MomentumIterativeFGSM(f, ϵ, μ, T::Int, args...; minval=0.0, maxval=1.0)
    x̂ = collect(args)
    F = eltype(first(x̂))
    α = F(ϵ/T); μ = F(μ)
    history = []
    g = [F(0.0) for i = 1:length(x̂)]
    for t = 1:T
        ∇x̂ = ∇xJ(f, x̂...)
        g = [(μ*gi .+ ∇x̂i/norm(∇x̂i,1)) for (∇x̂i,gi) in zip(∇x̂, g)]
        x̂ = [FGSM(x̂i,gi,α; minval=minval, maxval=maxval)
             for (x̂i,gi) in zip(x̂, g)]
        push!(history, x̂)
    end
    return history
end


text:
using Oceananigans.Solvers
using Oceananigans.Operators

"""
    implicit_free_surface_linear_operation!(result, x, arch, grid, bcs; args...)

Returns `L(ηⁿ)`, where `ηⁿ` is the free surface displacement at time step `n`
and `L` is the linear operator that arises
in an implicit time step for the free surface displacement `η`.

(See the docs section on implicit time stepping.)
"""
function implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)
    grid = L_ηⁿ⁺¹.grid
    arch = architecture(L_ηⁿ⁺¹)

    fill_halo_regions!(ηⁿ⁺¹, arch)

    event = launch!(arch, grid, :xy, _implicit_free_surface_linear_operation!,
                    L_ηⁿ⁺¹, grid,  ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt,
                    dependencies=Event(device(arch)))

    wait(device(arch), event)

    fill_halo_regions!(L_ηⁿ⁺¹, arch)

    return nothing
end

# Kernels that act on vertically integrated / surface quantities
@inline ∫ᶻ_Ax_∂x_ηᶠᶜᶜ(i, j, k, grid, ∫ᶻ_Axᶠᶜᶜ, η) = @inbounds ∫ᶻ_Axᶠᶜᶜ[i, j, k] * ∂xᶠᶜᵃ(i, j, k, grid, η)
@inline ∫ᶻ_Ay_∂y_ηᶜᶠᶜ(i, j, k, grid, ∫ᶻ_Ayᶜᶠᶜ, η) = @inbounds ∫ᶻ_Ayᶜᶠᶜ[i, j, k] * ∂yᶜᶠᵃ(i, j, k, grid, η)

"""
Compute the horizontal divergence of vertically-uniform quantity using
vertically-integrated face areas `∫ᶻ_Axᶠᶜᶜ` and `∫ᶻ_Ayᶜᶠᶜ`.
"""
@inline ∇h²ᶜᶜᵃ(i, j, k, grid, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, η::ReducedField{X, Y, Nothing}) where {X, Y} =
    1 / Azᶜᶜᵃ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, ∫ᶻ_Ax_∂x_ηᶠᶜᶜ, ∫ᶻ_Axᶠᶜᶜ, η) +
                                δyᵃᶜᵃ(i, j, k, grid, ∫ᶻ_Ay_∂y_ηᶜᶠᶜ, ∫ᶻ_Ayᶜᶠᶜ, η))

"""
    _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)

Returns the left hand side of the "implicit η equation"

```math
(∇ʰ⋅H∇ʰ - 1/gΔt²) ηⁿ⁺¹ = 1 / (gΔt) ∇ʰH U̅ˢᵗᵃʳ - 1 / (gΔt²) ηⁿ
----------------------
        ≡ L_ηⁿ⁺¹
```

which is written in a discrete finite volume form in which the equation
is arranged to ensure a symmtric form by multiplying by horizontal areas Az:

```
δⁱÂʷ∂ˣηⁿ⁺¹ + δʲÂˢ∂ʸηⁿ⁺¹ - 1/gΔt² Az ηⁿ⁺¹ = 1 / (gΔt) (δⁱÂʷu̅ˢᵗᵃʳ + δʲÂˢv̅ˢᵗᵃʳ) - 1 / gΔt² Az ηⁿ
```

where  ̂ indicates a vertical integral, and
       ̅ indicates a vertical average
"""
@kernel function _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)
    i, j = @index(Global, NTuple)

    @inbounds L_ηⁿ⁺¹[i, j, 1] = (  Azᶜᶜᵃ(i, j, 1, grid) * ∇h²ᶜᶜᵃ(i, j, 1, grid, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, ηⁿ⁺¹)
                                 - Azᶜᶜᵃ(i, j, 1, grid) * ηⁿ⁺¹[i, j, 1] / (g * Δt^2))
end


text:
@testset "Testing lss.jl" begin

    rough_kwargs = Dict(:atol => 1e-7, :rtol => 1e-7)

    # set up
    A = .95
    C = .05
    G = 1.
    mu_0 = [.75;]
    Sigma_0 = fill(0.000001, 1, 1)

    ss = LSS(A, C, G, mu_0)
    ss1 = LSS(A, C, G, mu_0, Sigma_0)

    vals = stationary_distributions(ss, max_iter=1000, tol=1e-9)

    @testset "test stationarity" begin
        vals = stationary_distributions(ss, max_iter=1000, tol=1e-9)
        ssmux, ssmuy, sssigx, sssigy = vals

        @test isapprox(ssmux, ssmuy; rough_kwargs...)
        @test isapprox(sssigx, sssigy; rough_kwargs...)
        @test isapprox(ssmux, [0.0]; rough_kwargs...)
        @test isapprox(sssigx, ss.C.^2 ./ (1 - ss.A .^2); rough_kwargs...)
    end

    @testset "test replicate" begin
        xval1, yval1 = replicate(ss, 100, 5000)
        xval2, yval2 = replicate(ss; t=100, num_reps=5000)
        xval3, yval3 = replicate(ss1; t=100, num_reps=5000)

        for (x, y)  in [(xval1, yval1), (xval2, yval2), (xval3, yval3)]
            @test isapprox(x , y; rough_kwargs...)
            @test abs(mean(x)) <= 0.05
        end

    end

    @testset "test convergence error" begin
        @test_throws ErrorException stationary_distributions(ss; max_iter=1, tol=eps())
    end

    @testset "test geometric_sums" begin
        β = 0.98
        xs = rand(10)
        for x in xs
            gsum_x, gsum_y = QuantEcon.geometric_sums(ss, β, [x])
            @test isapprox(gsum_x, ([x/(1-β *ss.A[1])]))
            @test isapprox(gsum_y, ([ss.G[1]*x/(1-β *ss.A[1])]))
        end
    end

    @testset "test constructors" begin
        # kwarg version
        other_ss = LSS(A, C, G; mu_0=[mu_0;])
        for nm in fieldnames(ss)
            @test getfield(ss, nm) == getfield(other_ss, nm)
        end
    end

end  # @testset


text:
module Filtering

include("../Types.jl")

export run

FILTER_RANGE = 1:12
FILTER_LENGTH = size(FILTER_RANGE)[1]
FILTER_STD = 0.35
HALF_FILTER_LENGTH = FILTER_LENGTH / 2
B = (FILTER_STD * HALF_FILTER_LENGTH) ^ 2

function generateCoefficients()
    coeff = []
    for i in FILTER_RANGE
        value = ℯ ^ (-0.5 * ((i - HALF_FILTER_LENGTH) / B))
        push!(coeff, value)
    end
    coeff
end

function run(input:: Channel)
    window = []
    filtersum = 0.0
    filtercoefficients = generateCoefficients()
    for i in FILTER_RANGE
        filtersum += filtercoefficients[i]
    end
    return function (output:: Channel)
        for dp in input
            push!(window, dp)
            if size(window)[1] == FILTER_LENGTH
                sum = 0.0
                for i in FILTER_RANGE
                    sum += window[i].magnitude * filtercoefficients[i]
                end
                new_dp = Types.Datapoint(window[Int(HALF_FILTER_LENGTH)].time, sum / filtersum)
                put!(output, new_dp)
                popfirst!(window)
            end
        end
    end
end

end

max_stars_repo_name: eliabiraschi/pedometer

max_issues_repo_name: eliabiraschi/pedometer

max_forks_repo_name: eliabiraschi/pedometer


text:
using TensorBoardLogger, Test
using TensorBoardLogger: preprocess, content
using TestImages
using ImageCore
@testset "TBText" begin
    data = Vector{Pair{String,Any}}()
    preprocess("test", TBText(Complex(1.0,3.0)), data)
    @test first(data) == ("test"=>TBText(Complex(1.0,3.0)))
end

@testset "TBHistogram" begin
    data = Vector{Pair{String,Any}}()
    arr = rand(ComplexF64, 5)
    preprocess("test", TBHistogram(arr), data)
    @test !(first(data) == ("test"=>TBHistogram(arr)))
    k,v = pop!(data)
    @test k == "test/im"
    @test content(v) == imag(arr) && typeof(v)<:TBHistogram
    k,v = pop!(data)
    @test k == "test/re"
    @test v.data == real(arr) && typeof(v)<:TBHistogram
end

@testset "TBVector" begin
    data = Vector{Pair{String,Any}}()
    arr = rand(ComplexF64, 5)
    preprocess("test", TBVector(arr), data)
    @test !(first(data) == ("test"=>TBVector(arr)))
    k,v = pop!(data)
    @test k == "test/im"
    @test content(v) == imag(arr) && typeof(v)<:TBVector
    k,v = pop!(data)
    @test k == "test/re"
    @test v.data == real(arr) && typeof(v)<:TBVector
end

@testset "TBImage" begin
    data = Vector{Pair{String,Any}}()
    testdata = [
    (rand(10).*100, L), (rand(10, 2), LN), (rand(2, 10), NL), (rand(10, 3), LC), (rand(3, 10), CL),
    (rand(10, 3, 2), LCN), (rand(3, 10, 2), CLN), (rand(2, 10, 3), NLC), (rand(2, 3, 10), NCL),
    (rand(4, 5), HW), (rand(5, 4), WH),
    (rand(5, 4, 2), WHN), (rand(4, 5, 2), HWN), (rand(2, 4, 5), NHW), (rand(2, 5, 4), NWH),
    (rand(5, 4, 3), WHC), (rand(4, 5, 3), HWC), (rand(3, 4, 5), CHW), (rand(3, 5, 4), CWH),
    (rand(5, 4, 3, 2), WHCN), (rand(4, 5, 3, 2), HWCN), (rand(2, 3, 4, 5), NCHW), (rand(2, 3, 5, 4), NCWH),
    (rand(3, 5, 4, 2), CWHN), (rand(3, 4, 5, 2), CHWN), (rand(2, 4, 5, 3), NHWC), (rand(2, 5, 4, 3), NWHC),
    ]
    for testtuple in testdata
        preprocess("randTBImage/$(last(testtuple))", TBImage(first(testtuple), last(testtuple)), data)
    end
    @test size(data, 1) == 45
end


text:
using ACE, LinearAlgebra
using ACE: evaluate

B = ACE.Utils.rpi_basis( species = [:X, :Al], N = 4,
                         pin = 0, pcut = 0)

Rs, Zs, z0 = ACE.Random.rand_nhd(10, B.pibasis.basis1p.J, :Al)
Zs[1] = 0

rcut = 3.
rcut_env = 6.
zcut_env = 3.

function fcut_env(r, z, rcut, zcut)
   return (r - rcut)^2 * (z - zcut)^2
end

function fcut_bond(r, rcut)
   return (r - rcut)^2
end

function envelope(Rs, rcut, rcut_env, zcut_env)
   fenv = fcut_bond(norm(Rs[1]), rcut)
   r̂ = Rs[1] / norm(Rs[1])
   o = Rs[1]/2
   for i = 2:length(Rs)
      R = Rs[i]
      z = dot(R - o, r̂)
      r = norm(R - z * r̂)
      fenv *= fcut_env(r, z, rcut_env, zcut_env)
   end
   return fenv
end

function eval_bond(B, Rs, Zs, z0)
   r̂ = Rs[1] / norm(Rs[1])
   o = Rs[1]/2
   Rr = map( r_ -> (r = r_ - o; o + r - 2*dot(r,r̂)*r̂), Rs )
   Rr[1] = Rs[1]
   return ( (evaluate(B, Rs, Zs, z0) + evaluate(B, Rr, Zs, z0))
            * envelope(Rs, rcut, rcut_env, zcut_env) )
end

b0 = eval_bond(B, Rs, Zs, z0)

# reversed configuration to test we have the right symmetry
Rs1 = [ r - Rs[1] for r in Rs ]
Rs1[1] = - Rs[1]
b1 = eval_bond(B, Rs1, Zs, z0)
b0 ≈ b1

B.pibasis.basis1p.J.
cutoff(B.pibasis.basis1p.J)
using Plots
rr = range(0, cutoff(B) + 1, length=100)
B1to5 = zeros(5, length(rr))
for (i, r) in enumerate(rr)
   B1to5[:, i] = evaluate(B.pibasis.basis1p.J, r)[1:5]
end
plot(rr, B1to5', ylims=[-2, 2])
vline!([B.pibasis.basis1p.J.rl, B.pibasis.basis1p.J.ru], lw=2, c=:black)


text:
# A simple grid world MDP
# All cells with reward are also terminal
# Environment created by: Zachary Sunberg

using CommonRLInterface
using StaticArrays
using Crayons
using Random: AbstractRNG

const RL = CommonRLInterface

# To avoid episodes of unbounded length, we put an arbitrary limit to the length of an
# episode. Because time is not captured in the state, this introduces a slight bias in
# the value function.
const EPISODE_LENGTH_BOUND = 200
const SIZE = SA[10, 10]
const REWARDS =     Dict(
		    SA[9,3] =>   1.0,
		    SA[8,8] =>   0.3,
		    SA[4,3] =>  -1.0,
		    SA[4,6] =>  -0.5)

mutable struct World <: AbstractEnv
  position ::SVector{2, Int}
  time :: Int
end

function World()
  return World(
    SA[rand(1:10), rand(1:10)],
    0)
end

RL.reset!(env::World, rng::AbstractRNG) = (env.position = SA[rand(rng, 1:SIZE[1]), rand(rng, 1:SIZE[2])])
RL.actions(env::World) = [SA[1,0], SA[-1,0], SA[0,1], SA[0,-1]]
RL.observe(env::World) = env.position

RL.terminated(env::World) =
  haskey(REWARDS, env.position) || env.time > EPISODE_LENGTH_BOUND

function RL.act!(env::World, a)
  # 40% chance of going in a random direction (=30% chance of going in a wrong direction)
  if rand() < 0.4
      a = rand(actions(env))
  end
  env.position = clamp.(env.position + a, SA[1,1], SIZE)
  env.time += 1
  return get(REWARDS, env.position, 0.0)
end

@provide RL.player(env::World) = 1 # An MDP is a one player game
@provide RL.players(env::World) = [1]
@provide RL.observations(env::World) = [SA[x, y] for x in 1:SIZE[1], y in 1:SIZE[2]]
@provide RL.clone(env::World) = World(env.position, env.time)
@provide RL.state(env::World) = env.position
@provide RL.setstate!(env::World, s) = (env.position = s)
@provide RL.valid_action_mask(env::World) = BitVector([1, 1, 1, 1])

# Additional functions needed by AlphaZero.jl that are not present in
# CommonRlInterface.jl. Here, we provide them by overriding some functions from
# GameInterface. An alternative would be to pass them as keyword arguments to
# CommonRLInterfaceWrapper.

function GI.render(env::World)
  for y in reverse(1:SIZE[2])
    for x in 1:SIZE[1]
      s = SA[x, y]
      r = get(REWARDS, s, 0.0)
      if env.position == s
        c = ("+",)
      elseif r > 0
        c = (crayon"green", "o")
      elseif r < 0
        c = (crayon"red", "o")
      else
        c = (crayon"dark_gray", ".")
      end
      print(c..., " ", crayon"reset")
    end
    println("")
  end
end

function GI.vectorize_state(env::World, state)
  v = zeros(Float32, SIZE[1], SIZE[2])
  v[state[1], state[2]] = 1
  return v
end

const action_names = ["r", "l", "u", "d"]

function GI.action_string(env::World, a)
  idx = findfirst(==(a), RL.actions(env))
  return isnothing(idx) ? "?" : action_names[idx]
end

function GI.parse_action(env::World, s)
  idx = findfirst(==(s), action_names)
  return isnothing(idx) ? nothing : RL.actions(env)[idx]
end

function GI.read_state(env::World)
  try
    s = split(readline())
    @assert length(s) == 2
    x = parse(Int, s[1])
    y = parse(Int, s[2])
    @assert 1 <= x <= SIZE[1]
    @assert 1 <= y <= SIZE[2]
    return SA[x, y]
  catch e
    return nothing
  end
end

GI.heuristic_value(::World) = 0.

GameSpec() = CommonRLInterfaceWrapper.Spec(World())


text:
using Pkg
pkg"activate ."
pkg"instantiate"
include(joinpath(@__DIR__, "src", "ReplitDatabase.jl"))
using .ReplitDatabase


text:
using DictWrappers
using FactCheck

facts("Failure") do
    for T in subtypes(DictWrapper)
        context(string(T)) do
            @fact_throws ArgumentError T(3)
            @fact_throws ArgumentError T(3.0)
            @fact_throws ArgumentError T((3, 3.0))
            @fact_throws ArgumentError T(())
            @fact_throws ArgumentError T(true)
        end
    end
end

facts("Simple composite") do
    test = 1 + im
    for T in subtypes(DictWrapper)
        context(string(T)) do
            wrapper = T(test)
            @fact typeof(wrapper) --> T
            @fact wrapper.wrapped --> exactly(test)

            fields = collect(keys(wrapper))
            @fact typeof(fields) --> Vector{Symbol}
            @fact length(fields) --> 2

            for field in fields
                @fact wrapper[field] --> exactly(getfield(test, field))
                value = wrapper[field]
                @fact value --> exactly(1)
            end

            for value in values(wrapper)
                @fact value --> exactly(1)
            end
        end
    end
end

facts("Complex composite") do
    test_part = 1 + im
    test = (test_part=>test_part)::Pair{Complex{Int64}, Complex{Int64}}

    context(string(RecursiveDictWrapper)) do
        wrapper = RecursiveDictWrapper(test)
        @fact typeof(wrapper) --> RecursiveDictWrapper
        @fact wrapper.wrapped --> exactly(test)

        fields = collect(keys(wrapper))
        @fact typeof(fields) --> Vector{Symbol}
        @fact length(fields) --> 2

        for field in fields
            @fact wrapper[field] --> RecursiveDictWrapper(getfield(test, field))
            @fact wrapper[field] --> RecursiveDictWrapper(test_part)
            value = wrapper[field]
            @fact value.wrapped --> exactly(test_part)

            subfields = collect(keys(value))
            @fact typeof(subfields) --> Vector{Symbol}
            @fact length(subfields) --> 2

            for subfield in subfields
                @fact value[subfield] --> getfield(test_part, subfield)
                subvalue = value[subfield]
                @fact subvalue --> exactly(1)
            end
        end

        @fact get(wrapper, :foo, 4) --> 4
        @fact (:first=>RecursiveDictWrapper(test_part)) in wrapper --> true
        @fact (:second=>RecursiveDictWrapper(test_part)) in wrapper --> true
        @fact (:first=>:wrong) in wrapper --> false
    end

    context(string(ShallowDictWrapper)) do
        wrapper = ShallowDictWrapper(test)
        @fact typeof(wrapper) --> ShallowDictWrapper
        @fact wrapper.wrapped --> exactly(test)

        fields = collect(keys(wrapper))
        @fact typeof(fields) --> Vector{Symbol}
        @fact length(fields) --> 2

        for field in fields
            @fact wrapper[field] --> getfield(test, field)
            @fact wrapper[field] --> test_part
        end

        @fact get(wrapper, :foo, 4) --> 4
        @fact (:first=>test_part) in wrapper --> true
        @fact (:second=>test_part) in wrapper --> true
        @fact (:first=>:wrong) in wrapper --> false
    end
end


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

using IterativeEigensolvers
using Test, LinearAlgebra, SparseArrays, Random

@testset "eigs" begin
    srand(1234)
    n = 10
    areal  = sprandn(n,n,0.4)
    breal  = sprandn(n,n,0.4)
    acmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4))
    bcmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4))

    testtol = 1e-6

    @testset for elty in (Float64, ComplexF64)
        if elty == ComplexF32 || elty == ComplexF64
            a = acmplx
            b = bcmplx
        else
            a = areal
            b = breal
        end
        a_evs = eigvals(Array(a))
        a     = convert(SparseMatrixCSC{elty}, a)
        asym  = copy(a') + a                  # symmetric indefinite
        apd   = a'*a                    # symmetric positive-definite

        b     = convert(SparseMatrixCSC{elty}, b)
        bsym  = copy(b') + b
        bpd   = b'*b

        (d,v) = eigs(a, nev=3)
        @test a*v[:,2] ≈ d[2]*v[:,2]
        @test norm(v) > testtol # eigenvectors cannot be null vectors
        (d,v) = eigs(a, LinearAlgebra.I, nev=3) # test eigs(A, B; kwargs...)
        @test a*v[:,2] ≈ d[2]*v[:,2]
        @test norm(v) > testtol # eigenvectors cannot be null vectors
        @test_logs (:warn,"Use symbols instead of strings for specifying which eigenvalues to compute") eigs(a, which="LM")
        @test_logs (:warn,"Adjusting ncv from 1 to 4") eigs(a, ncv=1, nev=2)
        @test_logs (:warn,"Adjusting nev from $n to $(n-2)") eigs(a, nev=n)
        # (d,v) = eigs(a, b, nev=3, tol=1e-8) # not handled yet
        # @test a*v[:,2] ≈ d[2]*b*v[:,2] atol=testtol
        # @test norm(v) > testtol # eigenvectors cannot be null vectors
        if elty <: LinearAlgebra.BlasComplex
            sr_ind = argmin(real.(a_evs))
            (d, v) = eigs(a, nev=1, which=:SR)
            @test d[1] ≈ a_evs[sr_ind]
            si_ind = argmin(imag.(a_evs))
            (d, v) = eigs(a, nev=1, which=:SI)
            @test d[1] ≈ a_evs[si_ind]
            lr_ind = argmax(real.(a_evs))
            (d, v) = eigs(a, nev=1, which=:LR)
            @test d[1] ≈ a_evs[lr_ind]
            li_ind = argmax(imag.(a_evs))
            (d, v) = eigs(a, nev=1, which=:LI)
            @test d[1] ≈ a_evs[li_ind]
        end

        (d,v) = eigs(asym, nev=3)
        @test asym*v[:,1] ≈ d[1]*v[:,1]
        @test eigs(asym; nev=1, sigma=d[3])[1][1] ≈ d[3]
        @test norm(v) > testtol # eigenvectors cannot be null vectors

        (d,v) = eigs(apd, nev=3)
        @test apd*v[:,3] ≈ d[3]*v[:,3]
        @test eigs(apd; nev=1, sigma=d[3])[1][1] ≈ d[3]

        (d,v) = eigs(apd, bpd, nev=3, tol=1e-8)
        @test apd*v[:,2] ≈ d[2]*bpd*v[:,2] atol=testtol
        @test norm(v) > testtol # eigenvectors cannot be null vectors

        @testset "(shift-and-)invert mode" begin
            (d,v) = eigs(apd, nev=3, sigma=0)
            @test apd*v[:,3] ≈ d[3]*v[:,3]
            @test norm(v) > testtol # eigenvectors cannot be null vectors

            (d,v) = eigs(apd, bpd, nev=3, sigma=0, tol=1e-8)
            @test apd*v[:,1] ≈ d[1]*bpd*v[:,1] atol=testtol
            @test norm(v) > testtol # eigenvectors cannot be null vectors
        end

        @testset "ArgumentErrors" begin
            @test_throws ArgumentError eigs(rand(elty,2,2))
            @test_throws ArgumentError eigs(a, nev=-1)
            @test_throws ArgumentError eigs(a, which=:Z)
            @test_throws ArgumentError eigs(a, which=:BE)
            @test_throws DimensionMismatch eigs(a, v0=zeros(elty,n+2))
            @test_throws ArgumentError eigs(a, v0=zeros(Int,n))
            if elty == Float64
                @test_throws ArgumentError eigs(a + copy(transpose(a)), which=:SI)
                @test_throws ArgumentError eigs(a + copy(transpose(a)), which=:LI)
                @test_throws ArgumentError eigs(a, sigma = rand(ComplexF32))
            end
        end
    end

    @testset "Symmetric generalized with singular B" begin
        srand(127)
        n = 10
        k = 3
        A = randn(n,n); A = A'A
        B = randn(n,k); B = B*B'
        @test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(eigvals(A, B)[1:k])
    end
end

# Problematic example from #6965A
let A6965 = [
    1.0   1.0   1.0   1.0   1.0   1.0   1.0  1.0
    -1.0   2.0   0.0   0.0   0.0   0.0   0.0  1.0
    -1.0   0.0   3.0   0.0   0.0   0.0   0.0  1.0
    -1.0   0.0   0.0   4.0   0.0   0.0   0.0  1.0
    -1.0   0.0   0.0   0.0   5.0   0.0   0.0  1.0
    -1.0   0.0   0.0   0.0   0.0   6.0   0.0  1.0
    -1.0   0.0   0.0   0.0   0.0   0.0   7.0  1.0
    -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  8.0
]
    d, = eigs(A6965,which=:SM,nev=2,ncv=4,tol=eps())
    @test d[1] ≈ 2.5346936860350002
    @test real(d[2]) ≈ 2.6159972444834976
    @test abs(imag(d[2])) ≈ 1.2917858749046127

    # Requires ARPACK 3.2 or a patched 3.1.5
    #T6965 = [ 0.9  0.05  0.05
    #          0.8  0.1   0.1
    #          0.7  0.1   0.2 ]
    #d,v,nconv = eigs(T6965,nev=1,which=:LM)
    # @test T6965*v ≈ d[1]*v atol=1e-6
end

# Example from Quantum Information Theory
import Base: size

mutable struct CPM{T<:LinearAlgebra.BlasFloat} <: AbstractMatrix{T} # completely positive map
    kraus::Array{T,3} # kraus operator representation
end
size(Phi::CPM) = (size(Phi.kraus,1)^2,size(Phi.kraus,3)^2)
LinearAlgebra.issymmetric(Phi::CPM) = false
LinearAlgebra.ishermitian(Phi::CPM) = false
function LinearAlgebra.mul!(rho2::StridedVector{T},Phi::CPM{T},rho::StridedVector{T}) where {T<:LinearAlgebra.BlasFloat}
    rho = reshape(rho,(size(Phi.kraus,3),size(Phi.kraus,3)))
    rho1 = zeros(T,(size(Phi.kraus,1),size(Phi.kraus,1)))
    for s = 1:size(Phi.kraus,2)
        As = view(Phi.kraus,:,s,:)
        rho1 += As*rho*As'
    end
    return copyto!(rho2,rho1)
end
LinearAlgebra.A_mul_B!(rho2::StridedVector{T},Phi::CPM{T},rho::StridedVector{T}) where {T<:LinearAlgebra.BlasFloat} = LinearAlgebra.mul!(rho2, Phi, rho)
# after the A_mul_B! deprecation, remove this A_mul_B! def

let
    # Generate random isometry
    (Q,R) = qr(randn(100,50))
    Q = reshape(Q,(50,2,50))
    # Construct trace-preserving completely positive map from this
    Phi = CPM(copy(Q))
    (d,v,nconv,numiter,numop,resid) = eigs(Phi,nev=1,which=:LM)
    # Properties: largest eigenvalue should be 1, largest eigenvector, when reshaped as matrix
    # should be a Hermitian positive definite matrix (up to an arbitrary phase)

    @test d[1] ≈ 1. # largest eigenvalue should be 1.
    v = reshape(v,(50,50)) # reshape to matrix
    v /= tr(v) # factor out arbitrary phase
    @test vecnorm(imag(v)) ≈ 0. # it should be real
    v = real(v)
    # @test vecnorm(v-v')/2 ≈ 0. # it should be Hermitian
    # Since this fails sometimes (numerical precision error),this test is commented out
    v = (v + v')/2
    @test isposdef(v)

    # Repeat with starting vector
    (d2,v2,nconv2,numiter2,numop2,resid2) = eigs(Phi,nev=1,which=:LM,v0=reshape(v,(2500,)))
    v2 = reshape(v2,(50,50))
    v2 /= tr(v2)
    @test numiter2 < numiter
    @test v ≈ v2

    # Adjust the tolerance a bit since matrices with repeated eigenvalues
    # can be very stressful to ARPACK and this may therefore fail with
    # info = 3 if the tolerance is too small
    @test eigs(sparse(1.0I, 50, 50), nev=10, tol = 5e-16)[1] ≈ fill(1., 10) #Issue 4246
end

@testset "real svds" begin
    A = sparse([1, 1, 2, 3, 4], [2, 1, 1, 3, 1], [2.0, -1.0, 6.1, 7.0, 1.5])
    S1 = svds(A, nsv = 2)
    S2 = svd(Array(A))

    ## singular values match:
    @test S1[1].S ≈ S2[2][1:2]
    @testset "singular vectors" begin
        ## 1st left singular vector
        s1_left = sign(S1[1].U[3,1]) * S1[1].U[:,1]
        s2_left = sign(S2[1][3,1]) * S2[1][:,1]
        @test s1_left ≈ s2_left

        ## 1st right singular vector
        s1_right = sign(S1[1].V[3,1]) * S1[1].V[:,1]
        s2_right = sign(S2[3][3,1]) * S2[3][:,1]
        @test s1_right ≈ s2_right
    end
    # Issue number 10329
    # Ensure singular values from svds are in
    # the correct order
    @testset "singular values ordered correctly" begin
        B = sparse(Diagonal([1.0, 2.0, 34.0, 5.0, 6.0]))
        S3 = svds(B, ritzvec=false, nsv=2)
        @test S3[1].S ≈ [34.0, 6.0]
        S4 = svds(B, nsv=2)
        @test S4[1].S ≈ [34.0, 6.0]
    end
    @testset "passing guess for Krylov vectors" begin
        S1 = svds(A, nsv = 2, v0=rand(eltype(A),size(A,2)))
        @test S1[1].S ≈ S2[2][1:2]
    end

    @test_throws ArgumentError svds(A,nsv=0)
    @test_throws ArgumentError svds(A,nsv=20)
    @test_throws DimensionMismatch svds(A,nsv=2,v0=rand(size(A,2)+1))

    @testset "Orthogonal vectors with repeated singular values $i times. Issue 16608" for i in 2:3
        rng = MersenneTwister(126) # Fragile to compute repeated values without blocking so we set the seed
        v0  = randn(rng, 20)
        d   = sort(rand(rng, 20), rev = true)
        for j in 2:i
            d[j] = d[1]
        end
        A = qr(randn(rng, 20, 20))[1]*Diagonal(d)*qr(randn(rng, 20, 20))[1]
        @testset "Number of singular values: $j" for j in 2:6
            # Default size of subspace
            F = svds(A, nsv = j, v0 = v0)
            @test F[1].U'F[1].U ≈ Matrix(I, j, j)
            @test F[1].V'F[1].V ≈ Matrix(I, j, j)
            @test F[1].S        ≈ d[1:j]
            for k in 3j:2:5j
                # Custom size of subspace
                F = svds(A, nsv = j, ncv = k, v0 = v0)
                @test F[1].U'F[1].U ≈ Matrix(I, j, j)
                @test F[1].V'F[1].V ≈ Matrix(I, j, j)
                @test F[1].S        ≈ d[1:j]
            end
        end
    end
end

@testset "complex svds" begin
    A = sparse([1, 1, 2, 3, 4], [2, 1, 1, 3, 1], exp.(im*[2.0:2:10;]), 5, 4)
    S1 = svds(A, nsv = 2)
    S2 = svd(Array(A))

    ## singular values match:
    @test S1[1].S ≈ S2[2][1:2]
    @testset "singular vectors" begin
        ## left singular vectors
        s1_left = abs.(S1[1].U[:,1:2])
        s2_left = abs.(S2[1][:,1:2])
        @test s1_left ≈ s2_left

        ## right singular vectors
        s1_right = abs.(S1[1].V[:,1:2])
        s2_right = abs.(S2[3][:,1:2])
        @test s1_right ≈ s2_right
    end
    @testset "passing guess for Krylov vectors" begin
        S1 = svds(A, nsv = 2, v0=rand(eltype(A),size(A,2)))
        @test S1[1].S ≈ S2[2][1:2]
    end

    @test_throws ArgumentError svds(A,nsv=0)
    @test_throws ArgumentError svds(A,nsv=20)
    @test_throws DimensionMismatch svds(A,nsv=2,v0=complex(rand(size(A,2)+1)))
end

@testset "promotion" begin
    eigs(rand(1:10, 10, 10))
    eigs(rand(1:10, 10, 10), rand(1:10, 10, 10) |> t -> t't)
    svds(rand(1:10, 10, 8))
    @test_throws MethodError eigs(big.(rand(1:10, 10, 10)))
    @test_throws MethodError eigs(big.(rand(1:10, 10, 10)), rand(1:10, 10, 10))
    @test_throws MethodError svds(big.(rand(1:10, 10, 8)))
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
#A julia script to load D4M if it isn't modularized.

include("../src/D4M.jl")


text:
module DataAndPlotting

using CSV
using HDF5
using Plots
# wave parameter
const ki = 2
const kj = 3

# Grid size
const ni = 30
const nj = 20

function index2coord(i , ni)
    @assert ni > 1
    @assert 1 <= i <= ni
    (i-1) / (ni-1)
end

# Set up standing wave
export initialize
function initialize()
    arr = zeros(Float64 , ni , nj)
    for j in 1:nj, i in 1:ni
        x = index2coord(i , ni)
        y = index2coord(j , nj)
        arr[i,j] = sin(2*π*ki*x) * sin(2*π*kj*y)
    end
    return arr
end

# Write data to file
export output
function output(arr::Array{Float64,2}, filename::String)
    h5write(filename, "data", arr)
end

# read file
export input
function input(filename::String)
    arr = h5read(filename, "data")
end

# Create a beautiful Plot
"not working"
export makeplot
function makeplot(arr::Array{Float64,2})
    x = range(0 , stop = 1 , length = ni)
    y = range(0 , stop = 1 , length = nj)
    p = plot(x, y, arr, st=[:surface, :contourf])
end

end # module


text:
using Test
using JuMP
using Cbc
using ConditionalJuMP
using Polyhedra
using LinearAlgebra

# In this example, we search for the point closest to the origin in the union
# of a list of convex polyhedra

# This implementation lets us add disjunctions and implications of the form
#   x in P
# where x is a vector and P is a Polyhedron
function ConditionalJuMP.Conditional(op::typeof(in), x::AbstractVector, P::HRepresentation)
    (&)([@?(h.a' * x <= h.β) for h in allhalfspaces(P)]...)
end

# A simple L1-norm objective we can minimize (since the Cbc
function l1norm_objective(x::AbstractVector{Variable})
    model = first(x).m
    y = @variable(model, [1:length(x)], basename="y", lowerbound=0)
    for i in 1:length(x)
        @constraint(model, y[i] >= x[i])
        @constraint(model, y[i] >= -x[i])
    end
    sum(y)
end

# P1 is a box from [-1.5, -1] to [-0.5, 1]
P1 = hrep(vcat(Matrix(I, 2, 2), -Matrix(I, 2, 2)), [-0.5, 1, 1.5, 1])

# P2 is a box from [1, -1] to [2, 1]
P2 = hrep(vcat(Matrix(I, 2, 2), -Matrix(I, 2, 2)), [2., 1, -1, 1])

m = Model(solver=CbcSolver(logLevel=0))
@variable(m, -5 <= x[1:2] <= 5)
@disjunction(m, x in P1, x in P2)
@objective(m, Min, l1norm_objective(x))

@test sum(m.colCat .== :Bin) == 1
solve(m)
x = getvalue(x)
@test x ≈ [-0.5, 0]


text:
using ViewerGL, LinearAlgebra
GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

# julia> tstarrangmnt(V,CV,FV,V); # pr-congrunc
faces2D =[
([-1.2798097711002834 -0.7175836704026446 -0.9273018142005801 -0.33575271353647373 -0.8888313802786274 -0.9234061838793065 -0.5515947470915735; -0.8529755793428437 -0.39350484584391116 -0.4216332551791834 -0.17487851279515618 -0.852196385854468 -0.4165123613451497 -0.17530866605859877], [[4, 7], [1, 5], [5, 4], [6, 1], [7, 2], [3, 6], [2, 3]]),
([0.4090297017807793 1.0911767094214293 1.1620015095819656 0.489477678628799 0.8016133264659618 0.7803006368816129; 0.6927213209798023 1.2289485323795242 1.6505941305087972 1.1716577255753056 1.0396095639207596 0.9571211278136456], [[1, 6], [3, 4], [4, 1], [2, 3], [6, 5], [5, 2]]),
([-1.0385711979250456 -0.7132409306523589 -1.0218142003120991 -1.0384325424416172 -0.7191489025549408; -0.41534838517272094 -0.3927141758263955 -0.1744935610128188 -0.4204698531446064 -0.17449353259972233], [[5, 3], [3, 1], [2, 5], [1, 4], [4, 2]]),
([-0.7258691288852619 -0.6448476135357281 -0.7529733231311402 -0.6448475748797735; -0.9414365540939522 -0.9488527179167867 -0.645323503925401 -0.6467259137115772], [[4, 2], [2, 1], [3, 4], [1, 3]]),
([-1.0121267711184847 -0.33480840100420095 -0.3348084695672869 -1.0121267493199877; -0.2597056379767005 0.21923078518196468 0.6981670688632455 0.2192306971122646], [[2, 1], [4, 3], [3, 2], [1, 4]]),
([0.6448476063306642 1.3221658778195489 1.322165901959854 0.6448475232354451; -0.1463457735868241 0.3325905535984814 0.8115269773474794 0.3325905612647395], [[4, 1], [2, 3], [3, 4], [1, 2]]),
([-0.06604591075614946 0.18807124440577458 0.1905165467536977 -0.04021340022048541 0.14217000611740668; -0.04713960699351189 0.2068955875876711 0.209505626479912 -0.02326566448216332 0.16361675908761078], [[4, 2], [1, 4], [2, 3], [3, 5], [5, 1]]),
([-0.6697247396897046 -0.2086932468973245 -0.5052490695995919 0.006994807111632216 -0.24632547950176553 -0.24927887576058852 0.15179576190233096; 0.644847403065945 0.7568013675147519 1.322165656894694 0.7258689440961424 0.6448474196145497 0.6502453745761758 1.3221657898108872], [[5, 1], [3, 7], [7, 4], [1, 3], [6, 5], [2, 6], [4, 2]]),
([-0.6503139392786357 -0.29398729456104555 -0.37586972830149784 -0.5732109629805938 -0.41272353233912057 -0.5233786186832611 -0.46562652421258094; -0.7396197100048467 -0.3972361696037834 -0.45399913155911314 -0.6186936367403125 -0.5834586434704805 -0.6758077398003293 -0.5948847636314626], [[1, 7], [7, 6], [2, 3], [5, 2], [4, 1], [6, 5], [3, 4]]),
([0.6467332598875593 0.7851007085073018 0.7851007074696446 0.6776070311248275; 0.018049561024192355 0.09608628388193051 0.2633898315127631 0.21470523454981943], [[4, 1], [3, 4], [1, 2], [2, 3]]),
([-0.779895602742659 -0.4553408517346448 -0.652509317340201 -0.7650426585272763 -0.6534665149433275 -0.6579726462257522 -0.4560220120135956 -0.4411690677982131; -0.1673539149934392 -0.046897891072995655 0.17587361979499583 0.13889993356702787 0.1708238266062496 0.17407860148447787 -0.06094278034815384 0.24531106821231316], [[3, 5], [5, 6], [1, 7], [6, 4], [8, 3], [4, 1], [7, 2], [2, 8]]),
([0.7622970369442709 0.7692584135831895 0.76553099855101; 0.31581372865520774 0.3201890192230694 0.32095917543894364], [[1, 3], [2, 1], [3, 2]]),
([-1.0541060033025456 -0.7283956506338671 -1.0546135163471562 -0.7416434996459378 -1.0535823727283933; 0.5846207610275508 0.583084917238596 0.6922503610772871 0.7208776421538298 0.5791744063734227], [[2, 5], [1, 3], [4, 2], [5, 1], [3, 4]]),
([0.46088661668800884 0.6467332598875607 0.4608866167915211 0.6724489914908843; 0.6695507645610579 0.6695507646203934 0.9937648018995732 0.9937648018377638], [[4, 2], [2, 1], [3, 4], [1, 3]]),
([-0.36557656064828625 0.5922962739322862 0.11335982892788275 0.11335982101674907; -0.8518116882414981 -0.17449331351230066 -0.17449334847177042 -0.8518116436467215], [[3, 2], [4, 1], [1, 3], [2, 4]]),
([-0.3899181623603573 -0.06555604191791929 -0.3896220982289381 -0.15929960118731126 -0.07535116145235445 -0.1969108414411048; 0.9839649649611354 0.9937645583723279 0.9937645541302548 0.9642674642570118 0.6695505212320305 0.6695505196407963], [[6, 4], [4, 1], [5, 6], [3, 2], [2, 5], [1, 3]]),
([0.2696522928626458 0.5018086808859741 0.47513192706916185 0.2732636209749537; -0.6421317134295561 -0.514805620005238 -0.26747630378962917 -0.3925646287592917], [[1, 4], [2, 1], [4, 3], [3, 2]]),
([-0.6922503610772871 -0.5846207610275508 -0.6913581132960885; 0.9839649649611354 0.9937645871879821 0.9937645541302549], [[1, 2], [2, 3], [3, 1]]),
([-0.0003558289156989935 0.7475543240956147 0.7475543193661965 -0.00035575654402898125; -0.3986377043306246 0.13021464464356017 0.6590669546828719 0.1302146396995053], [[3, 2], [1, 4], [4, 3], [2, 1]]),
([-0.5634030165076495 0.18450710897255157 -0.5634030323319841 0.1845070931482169; -0.8054012465702503 -0.8054012623945849 -0.057491114354159815 -0.05749113017849419], [[4, 3], [1, 2], [2, 4], [3, 1]]),
([0.057491069890257734 0.8054012021063481 0.8054012179306865 0.057491085714595935; -0.3577337617253023 0.17111858272300895 0.699970876589569 0.17111853214125755], [[4, 1], [2, 3], [3, 4], [1, 2]]),
([-0.8054012623945846 -0.057491130178494144 -0.8054012179306864 -0.05749108571459594; -0.7475543709919092 -0.7475543265280078 0.0003557049183168017 0.0003557493822183141], [[2, 4], [3, 1], [1, 2], [4, 3]]),
([-0.9453779886878261 -0.047691586676358796 -0.9347162589906084 -0.6706800769589929 -0.6181290860595076 -0.6202522863734578 -1.0845096880603842 -0.46152119777774997 -0.9453833058121955; -0.48963784662390314 0.0008336133858533825 -0.4701241938398577 0.00704930012507864 -0.3787325971054617 -0.28454624337107204 -0.7408350235914223 -0.7470507103306474 -0.48940197600963187], [[8, 2], [9, 7], [4, 3], [7, 8], [2, 4], [5, 1], [3, 6], [6, 5], [1, 9]]),
([-0.4768341894077711 0.0003557844517982159 0.0003557120801280599 -0.7475543685595155 -0.7475543638300972 -0.47941960720940857 -0.4961126373643564; -0.7314398673985791 -0.4747709995961227 0.15573526438373286 -0.8770536599771404 -0.24654733005280857 -0.7278644368281862 -0.7418092751419646], [[5, 4], [2, 3], [3, 5], [4, 7], [1, 2], [7, 6], [6, 1]]),
([-0.32426714810061874 -0.3178339849732556 -0.31651935101416007; 0.47683418940777095 0.47941960720940846 0.49611263736435623], [[3, 1], [2, 3], [1, 2]]),
([-0.48963784662390314 -0.4701241938398577 -0.48940197600963187 -0.28454624337107204 -0.3787325971054617; 0.049730280157139754 0.06587213039119212 0.050015420782415755 0.07643840590530138 -0.03742213531979252], [[2, 3], [1, 5], [4, 2], [5, 4], [3, 1]]),
([0.6695507643553293 0.9837550033206118 0.9937648017030808 0.6695507644234329; 0.6733097772450742 0.9871959122403988 0.9978523852459957 0.6781379408931124], [[3, 2], [1, 4], [2, 1], [4, 3]]),
([0.9837550033206085 0.9937648016938412 0.9937648017030773; 0.6376220802367375 0.6432822564072982 0.6553240139833155], [[1, 3], [2, 1], [3, 2]]),
([-0.818745341363344 -0.5007131023573861 -0.5215884382582039 -0.49437881059673994 -0.49472235965758204 -0.8076180014889325; -1.2645868510036822 -1.0460582901346591 -0.8321240422303853 -1.0425221183871403 -1.050026387139559 -1.0215261503271076], [[2, 4], [3, 2], [6, 3], [5, 1], [1, 6], [4, 5]]),
([0.6787799844049875 0.992143482099084 0.9929193167327109 0.6795526588743764; 0.6695507643553293 0.9828495290877188 0.9837550033206118 0.6695507643674667], [[3, 2], [2, 4], [4, 1], [1, 3]]),
([0.46088661502747014 0.4989025303730048 0.4608866155771692; 0.09562318959388295 0.09562318978859051 0.20294985482814365], [[3, 2], [2, 1], [1, 3]]),
([0.6695507645610579 0.9937648018995732 0.9937648017030811 0.6695507644234331; 0.19436405547083296 0.42883649062481927 0.6426672510025468 0.3441542736299178], [[4, 1], [1, 2], [2, 3], [3, 4]])
];

for (V,EV) in faces2D
	VV = [[k] for k=1:size(V,2)]
	model = (V, [VV,EV])::Lar.LARmodel
	println(EV,“\n”)
	GL.VIEW(GL.numbering(.1)(model, GL.COLORS[1], 0.1));
end


text:
# if VERSION >=v"1.3.0"
# 	(G::Group)(gens::AbstractVector{<:Perm}) = PermGroup(gens)
# end

# embeddings

function emb!(res::Perm, g::AbstractPerm, V::AbstractVector{<:Integer})
    @assert degree(res) >= degree(g)
	for (idx, v) in pairs(V)
		res[v] = idx^g
	end
    return res
end

emb(g::AbstractPerm, n::Integer) = emb!(Perm(n), g, eachindex(g))

function (G::PermGroup)(p::AbstractPerm)
	parent(p) === G && return p
    g = degree(p) == degree(G) ? p : emb(p, degree(G))
	q = perm(g)
	q in G && return Permutation(q, G)
	throw("Cannot coerce: $p is not an element of $G")
end

@doc doc"""
    base(G::PermGroup)
Get the base of the group `G`.
"""
base(G::AbstractPermutationGroup) = base(StabilizerChain(G))

@doc doc"""
    sgs(G::PermGroup)
Compute strong generating set for group `G`.

If a stabilizer chain for `G` has been already stored in `G`, `sgs` will be
just retrieved from it.
"""
sgs(G::AbstractPermutationGroup) = sgs(StabilizerChain(G))

@doc doc"""
    in(g::perm, G::PermGroup)
Membership test for permutation group `G` by `sift`ing `g` through `StabilizerChain(G)`.
"""
function Base.in(g::AbstractPerm, G::AbstractPermutationGroup)
    g = degree(g) < degree(G) ? emb(g, degree(G)) : g
    return g in StabilizerChain(G)
end

@doc doc"""
	transversals(G::PermGroup)
Return the transversals (as a Vector) of a permutation group `G`.
"""
transversals(G::AbstractPermutationGroup) = transversals(StabilizerChain(G))

@doc doc"""
    perm_by_baseimages(G::PermGroup, base_images::AbstractVector{<:Integer})
Return the unique permutation in `G` determined by to `base_images`
(with respect to `base(G)`).
"""
function perm_by_baseimages(
    G::AbstractPermutationGroup,
    baseimages::AbstractVector{<:Integer},
	inverse = true,
	tmp_word = Word(eltype(eltype(G))[])
)
    @boundscheck length(baseimages) == length(base(G))
    trans = transversals(G)
    res = one(G)

	for (i, schr) in enumerate(trans)
		bi = baseimages[i]
	    bi == first(schr) && continue
	    word = append!(tmp_word, schr.gens_inv, schr.orb, bi, schr.op)
	    res = tmp_word(schr.gens_inv, res)
		empty!(tmp_word)
	end

    return inverse ? inv(res) : res
end

###

struct BaseImagesIter{I, T}
	images::Vector{I}
	transversals::T
end

Base.eltype(::Type{<:BaseImagesIter{I}}) where {I} = Vector{I}
Base.length(bitr::BaseImagesIter) = prod(length, bitr.transversals)

base_images(G::PermGroup) = BaseImagesIter(copy(base(G)), transversals(G))

@inline Base.iterate(bitr::BaseImagesIter) = bitr.images, (count=1, total=length(bitr))

@inline function Base.iterate(bitr::BaseImagesIter, state)
	state.count == state.total && return nothing
	img = next!(bitr.images, bitr.transversals)
	return img, (count=state.count+1, total=state.total)
end

@inline function next!(baseimages::AbstractVector{<:Integer}, transversals)
	for position in length(baseimages):-1:1
		t = transversals[position]
		if islast(t, baseimages[position])
			# @debug "last point in orbit: spilling at" position collect(t), baseimages[position]
			baseimages[position] = first(t)
		else
			# @debug "next point in orbit: incrementing at" position collect(t), baseimages[position]
			baseimages[position] = _next(t, baseimages[position])
			break
		end
	end
	return baseimages
end

function _next(itr, elt)
	x, state = iterate(itr)
	while x != elt
		x, state = iterate(itr, state)
	end
	return first(iterate(itr, state))
end


text:
module Zygote

using LinearAlgebra

# This flag enables Zygote to grab extra type inference information during
# compiles. When control flow is present, this can give gradient code a
# performance boost.

# HOWEVER, this is not Jameson-approved, nor well supported by the compiler, and
# has several caveats. Recursion will cause inference to stack overflow.
# Gradient redefinitions may result in ugly type errors. And Jameson *will* know.
const usetyped = get(ENV, "ZYGOTE_TYPED", false) == "true"

using IRTools
using MacroTools, Requires
using MacroTools: @forward

export Params, gradient, derivative, forward, @code_grad

include("tools/idset.jl")
include("tools/ir.jl")
include("tools/reflection.jl")
include("tools/fillarray.jl")

include("compiler/reverse.jl")
include("compiler/emit.jl")
include("compiler/interface.jl")

include("forward/Forward.jl")

include("lib/grad.jl")
include("lib/lib.jl")
include("lib/real.jl")
include("lib/complex.jl")
include("lib/base.jl")
include("lib/array.jl")
include("lib/nnlib.jl")
include("lib/broadcast.jl")

# we need to define this late, so that the genfuncs see lib.jl
include("compiler/interface2.jl")
usetyped || include("precompile.jl")

@init @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("flux.jl")

# helps to work around 265-y issues
function refresh()
  include(joinpath(@__DIR__, "compiler/interface2.jl"))
  include(joinpath(@__DIR__, "precompile.jl"))
  return
end

end # module


text:
#---

using ACE1
using Printf, Test, LinearAlgebra, JuLIP, JuLIP.Testing, Random
using JuLIP: evaluate, evaluate_d, evaluate_ed
using ACE1: combine
using BenchmarkTools

#---

maxdeg = 13
N = 4
species = [:Ti, :Al]
r0 = 1.0
rcut = 3.0
trans = PolyTransform(1, r0)
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
D = ACE1.SparsePSHDegree()
P1 = ACE1.BasicPSH1pBasis(Pr; species = species)
basis = ACE1.PIBasis(P1, N, D, maxdeg)

#---

c = randcoeffs(basis)
p = 0.05
csp = copy(c)
for i = 1:length(csp)
   if rand() > p
      csp[i] = 0
   end
end

#---

Vdag = combine(basis, csp)
V = standardevaluator(Vdag)
Vdag_sp = ACE1.deletezeros(Vdag)
V_sp = standardevaluator(Vdag_sp)
# Vdag_sp = graphevaluator(V_sp)

#---

Nat = 3
Rs, Zs, z0 = ACE1.rand_nhd(Nat, Pr, species)

# evaluate_d(Vdag_sp, Rs, Zs, z0)

val0 = evaluate(V, Rs, Zs, z0)
dV0 = evaluate_d(V, Rs, Zs, z0)
for (V, str) in ( (Vdag, "Vdag"), (V_sp, "V_sp"), (Vdag_sp, "Vdag_sp") )
   @info("Error for $str:")
   val1 = evaluate(V, Rs, Zs, z0)
   @show abs(val0 - val1)
   dV1 = evaluate_d(V, Rs, Zs, z0)
   @show norm(dV0 - dV1)
end

#---

for (pot, name) in ( (V, "V"), (Vdag, "Vdag"),
                   (V_sp, "V_sp"), (Vdag_sp, "Vdag_sp") )
   tmp = ACE1.alloc_temp(pot, Nat)
   @info("Energy Timing for $name")
   @btime ACE1.evaluate!($tmp, $pot, $Rs, $Zs, $z0)
end

#---

graphevaluator(Vdag_sp)

for (pot, name) in ( (V, "V"), (Vdag, "Vdag"),
                   (V_sp, "V_sp"), (Vdag_sp, "Vdag_sp") )
   tmpd = ACE1.alloc_temp_d(pot, Nat)
   @info("Force Timing for $name")
   @btime ACE1.evaluate_d!($(tmpd.dV), $tmpd, $pot, $Rs, $Zs, $z0)
end


text:
"""
`XFloats.jl` -- fast 32-bit and 16-bit floats with extended accuracy.

Two types are exported:
- `XFloat32`. 32-bit floating point type with eXtended computational accuracy.
- `XFloat16`. 16-bit floating point type with eXtended computational accuracy.
"""
module XFloats

export XFloat32, XFloat16

using Printf
using Random
using LinearAlgebra
using SpecialFunctions

"""
    XFloat32 (a float type)

A more accurate alternative to Float32 that works as a Float64.
"""
primitive type XFloat32 <: AbstractFloat 64 end
XFloat32(x::XFloat32) = x

"""
    XFloat16 (a float type)

A more accurate alternative to Float16 that works as a Float32.
"""
primitive type XFloat16 <: AbstractFloat 32 end
XFloat16(x::XFloat16) = x

include("type/construct.jl")
include("type/promote_convert.jl")
include("type/representations.jl")

include("imports.jl")
include("dispatch/specialize.jl")
include("dispatch/linearalgebra.jl")

end # module XFloats


text:
using FileIO
using FilePathsBase
using Test

# Both FileIO and FilePathsBase export filename, but we only want the FileIO definition.
using FileIO: filename

struct MimeSaveTestType
end

@testset "FileIO" begin
    include("query.jl")
    include("loadsave.jl")
    include("error_handling.jl")
    include("test_mimesave.jl")
end


text:
using TensorOperations
using CUDA

let
    n = 1 << 12
    # test a common operation when contracting two tensors
    a = CUDA.rand(ComplexF32, n, 2, n)
    b = CUDA.rand(ComplexF32, n, 2, n)

    a_indices = Tuple([-1, -2, 1])
    b_indices = Tuple([1, -3, -4])
    c_indices = Tuple(symdiff(a_indices, b_indices))

    CUDA.@profile c = tensorcontract(a, a_indices, b,b_indices, c_indices)

    @show typeof(c)
    @show size(c)
end


text:
using Test
using LinearAlgebra
using QuantumControlBase.TestUtils
using QuantumControlBase: AbstractControlObjective, WeightedObjective

@testset "Sparse objective adjoint" begin

    obj = dummy_control_problem().objectives[1]
    adj = adjoint(obj)

    @test norm(adj.initial_state - obj.initial_state') ≈ 0
    @test norm(adj.target_state - obj.target_state') ≈ 0
    @test norm(adj.generator[1] - obj.generator[1]') ≈ 0
    @test norm(adj.generator[2][1] - obj.generator[2][1]') ≈ 0
    @test adj.generator[2][2] == obj.generator[2][2]

end

@testset "Dense objective adjoint" begin

    obj = dummy_control_problem(sparsity=1.0).objectives[1]
    adj = adjoint(obj)

    @test norm(adj.initial_state - obj.initial_state') ≈ 0
    @test norm(adj.target_state - obj.target_state') ≈ 0
    @test norm(adj.generator[1] - obj.generator[1]') ≈ 0
    @test norm(adj.generator[1] - obj.generator[1]) ≈ 0
    @test norm(adj.generator[2][1] - obj.generator[2][1]') ≈ 0
    @test adj.generator[2][2] == obj.generator[2][2]

end

@testset "Non-Hermitian objective adjoint" begin

    obj = dummy_control_problem(sparsity=1.0, hermitian=false).objectives[1]
    adj = adjoint(obj)

    @test norm(adj.initial_state - obj.initial_state') ≈ 0
    @test norm(adj.target_state - obj.target_state') ≈ 0
    @test norm(adj.generator[1] - obj.generator[1]') ≈ 0
    @test norm(adj.generator[1] - obj.generator[1]) > 0
    @test norm(adj.generator[2][1] - obj.generator[2][1]') ≈ 0
    @test norm(adj.generator[2][1] - obj.generator[2][1]) > 0
    @test adj.generator[2][2] == obj.generator[2][2]

end

@testset "weighted objective adjoint" begin

    obj0 = dummy_control_problem().objectives[1]
    obj = WeightedObjective(
        initial_state=obj0.initial_state,
        generator=obj0.generator,
        target_state=obj0.target_state,
        weight=0.2
    )
    adj = adjoint(obj)

    @test adj.weight == obj.weight

end

@testset "custom objective adjoint" begin

    obj0 = dummy_control_problem(hermitian=false).objectives[1]

    Base.@kwdef struct _CustomGateObjective <: AbstractControlObjective
        initial_state
        generator
        gate
        weight
        coeff
    end

    obj = _CustomGateObjective(obj0.initial_state, obj0.generator, "CNOT", 0.5, 1im)
    adj = adjoint(obj)

    @test norm(adj.initial_state - obj.initial_state') ≈ 0
    @test norm(adj.generator[1] - obj.generator[1]') ≈ 0
    @test norm(adj.generator[1] - obj.generator[1]) > 0
    @test adj.gate == "CNOT"
    @test adj.weight == 0.5
    @test adj.coeff == -1im

end