text:
using JuMP, Ipopt, Optim, LineSearches, LinearAlgebra, SparseArrays, Distributions, DataFrames, CSV, HTTP

function dense_hessian(hessian_sparsity, V, n)
    I = [i for (i,j) in hessian_sparsity]
    J = [j for (i,j) in hessian_sparsity]
    raw = sparse(I, J, V, n, n)
    return Matrix(raw + raw' - sparse(diagm(0=>diag(raw))))
end

function wrapper()
    # Let's read in the data from PS8
    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS8-factor/nlsy.csv"
    df = CSV.read(HTTP.get(url).body)
    X = [df.black df.hispanic df.female df.schoolt df.gradHS df.grad4yr ones(size(df,1),1)]
    y = df.logwage

    # first let's do unconstrained optimization
    function reg_mle(θ, X, y)
        # first K elements are the coefficients of the outcome equation
        β = θ[1:end-1]
        # last element is the variance (stdev)
        σ = θ[end]
        # now build the likelihood
        loglike = -sum(-.5 .* ( log(2*π) .+ log(σ^2) .+ ( (y .- X*β)./σ ).^2 ) )
        # more intuitive way? (but JuMP can't use pdf's from Distributions.jl)
        #loglike = -sum( log(1 ./ sqrt(σ^2)) .+ logpdf.(Normal(0,1),(y .- X*β)./sqrt(σ^2)) )
        return loglike
    end

    # run the optimizer for MLE
    svals = vcat(X\y,.5);
    td = TwiceDifferentiable(th -> reg_mle(th, X, y), svals; autodiff = :forward)
    θ̂_optim_ad = optimize(td, svals, Newton(linesearch = BackTracking()), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=1))
    θ̂_mle_optim_ad = θ̂_optim_ad.minimizer
    loglikeval = θ̂_optim_ad.minimum
    # evaluate the Hessian at the estimates
    H  = Optim.hessian!(td, θ̂_mle_optim_ad)
    θ̂_mle_optim_ad_se = sqrt.(diag(inv(H)))
    # store results in a data frame
    results = DataFrame(coef_mle = vcat(vec(θ̂_mle_optim_ad),-loglikeval), se_mle = vcat(vec(θ̂_mle_optim_ad_se),missing), coef_ols = vcat(X\y,missing,missing) )

    # now use JuMP (objective function + optimization call is all in one function here)
    function jump_mle(θ₀, X, y)
        # define the model
        model = Model(Ipopt.Optimizer)
        set_silent(model)
        @variable(model, β[j=1:size(X,2)], start = θ₀[j])
        @variable(model, σ, start = θ₀[end])
        @NLobjective(model, Max, sum(-.5 * ( log(2*π) + log(σ^2) + ((y[i] - sum(X[i,j]*β[j] for j in 1:size(X,2)) )/σ)^2 ) for i in 1:size(X,1) ) )
        # optimize the model
        JuMP.optimize!(model)
        # return parameter estimates
        coef_jump = vcat(JuMP.value.(β), JuMP.value(σ), JuMP.objective_value(model) )
        # return Hessian for SEs
        values = coef_jump[1:end-1]
        MOI = JuMP.MathOptInterface
        d = JuMP.NLPEvaluator(model)
        MOI.initialize(d, [:Hess])
        hessian_sparsity = MOI.hessian_lagrangian_structure(d)
        V = zeros(length(hessian_sparsity))
        MOI.eval_hessian_lagrangian(d, V, values, 1.0, Float64[])
        H = dense_hessian(hessian_sparsity, V, length(values))
        se_jump = sqrt.(diag(inv(-H)))
        return coef_jump, se_jump
    end
    jump_coefs,jump_se = jump_mle(svals, X, y)
    results.coef_jump = jump_coefs
    results.se_jump = vcat(jump_se,missing)

    # now let's use JuMP to estimate a constrained optimization
    function jump_cns_mle(θ₀, X, y)
        # define the model
        model = Model(Ipopt.Optimizer)
        set_silent(model)
        @variable(model, β[j=1:size(X,2)], start = θ₀[j])
        @variable(model, σ, start = θ₀[end])
        @constraint(model, β[2] == .16)
        @NLobjective(model, Max, sum(-.5 * ( log(2*π) + log(σ^2) + ((y[i] - sum(X[i,j]*β[j] for j in 1:size(X,2)) )/σ)^2 ) for i in 1:size(X,1) ) )
        # optimize the model
        JuMP.optimize!(model)
        # return parameter estimates
        coef_jump = vcat(JuMP.value.(β), JuMP.value(σ), JuMP.objective_value(model) )
        # return Hessian for SEs
        values = coef_jump[1:end-1]
        MOI = JuMP.MathOptInterface
        d = JuMP.NLPEvaluator(model)
        MOI.initialize(d, [:Hess])
        hessian_sparsity = MOI.hessian_lagrangian_structure(d)
        V = zeros(length(hessian_sparsity))
        MOI.eval_hessian_lagrangian(d, V, values, 1.0, Float64[])
        H = dense_hessian(hessian_sparsity, V, length(values))
        se_jump = sqrt.(diag(inv(-H)))
        return coef_jump, se_jump
    end
    jump_cns_coefs,jump_cns_se = jump_cns_mle(svals, X, y)
    results.coef_jump_cns = jump_cns_coefs
    results.se_jump_cns = vcat(jump_cns_se,missing)
    @show results

    # can we use optim for the same constrained optimization?
    function cns_reg_mle(θ, cns_mat, X, y)
        # first K elements are the coefficients of the outcome equation
        β = θ[1:end-1]
        # last element is the variance (stdev)
        σ = θ[end]
        # impose constraint
        for r in size(cns_mat,1)
            insert!(β,convert(Int64,cns_mat[r,1]),cns_mat[r,2])
        end
        # now build the likelihood
        loglike = -sum(-.5 .* ( log(2*π) .+ log(σ^2) .+ ( (y .- X*β)./σ ).^2 ) )
        # more intuitive way? (but JuMP can't use pdf's from Distributions.jl)
        #loglike = -sum( log(1 ./ sqrt(σ^2)) .+ logpdf.(Normal(0,1),(y .- X*β)./sqrt(σ^2)) )
        return loglike
    end

    # run the optimizer for MLE
    # first, set up constraints
    cns_mat = hcat(2,.16)
    svals = vcat(X\y,.5)
    # constraints are treated as data, so take them out of starting values (they get added back in inside the obj function)
    deleteat!(svals, convert(Int64,cns_mat[1,1]))
    td = TwiceDifferentiable(th -> cns_reg_mle(th, cns_mat, X, y), svals; autodiff = :forward)
    θ̂_optim_ad = optimize(td, svals, Newton(linesearch = BackTracking()), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=1))
    θ̂_mle_optim_ad = θ̂_optim_ad.minimizer
    loglikeval = θ̂_optim_ad.minimum
    # evaluate the Hessian at the estimates
    H  = Optim.hessian!(td, θ̂_mle_optim_ad)
    θ̂_mle_optim_ad_se = sqrt.(diag(inv(H)))
    # add back in constraint in both estimates and SEs
    insert!(θ̂_mle_optim_ad   ,convert(Int64,cns_mat[1,1]),cns_mat[1,2])
    insert!(θ̂_mle_optim_ad_se,convert(Int64,cns_mat[1,1]),0)
    # store results in a data frame
    results.coef_optim_cns = vcat(θ̂_mle_optim_ad,-loglikeval)
    results.se_optim_cns   = vcat(θ̂_mle_optim_ad_se,missing)
    @show results

    # now let's use JuMP to estimate another constrained optimization
    function jump_cns2_mle(θ₀, X, y)
        # define the model
        model = Model(Ipopt.Optimizer)
        set_silent(model)
        @variable(model, β[j=1:size(X,2)], start = θ₀[j])
        @variable(model, σ, start = θ₀[end])
        @constraint(model, β[2] == .16)
        @constraint(model, β[4] == 1+2*β[3])
        @NLobjective(model, Max, sum(-.5 * ( log(2*π) + log(σ^2) + ((y[i] - sum(X[i,j]*β[j] for j in 1:size(X,2)) )/σ)^2 ) for i in 1:size(X,1) ) )
        # optimize the model
        JuMP.optimize!(model)
        # return parameter estimates
        coef_jump = vcat(JuMP.value.(β), JuMP.value(σ), JuMP.objective_value(model) )
        # return Hessian for SEs
        values = coef_jump[1:end-1]
        MOI = JuMP.MathOptInterface
        d = JuMP.NLPEvaluator(model)
        MOI.initialize(d, [:Hess])
        hessian_sparsity = MOI.hessian_lagrangian_structure(d)
        V = zeros(length(hessian_sparsity))
        MOI.eval_hessian_lagrangian(d, V, values, 1.0, Float64[])
        H = dense_hessian(hessian_sparsity, V, length(values))
        se_jump = sqrt.(diag(inv(-H)))
        return coef_jump, se_jump
    end
    jump_cns2_coefs,jump_cns2_se = jump_cns2_mle(svals, X, y)
    results.coef_jump_cns2 = jump_cns2_coefs
    results.se_jump_cns2 = vcat(jump_cns2_se,missing)
    @show results

    # can we use optim for the same constrained optimization?
    function cns2_reg_mle(θ, cns_mat, X, y)
        # first K elements are the coefficients of the outcome equation
        β = θ[1:end-1]
        # last element is the variance (stdev)
        σ = θ[end]
        # impose constraints
        for r in 1:size(cns_mat,1)
            idx1 = convert(Int64,cns_mat[r,1])
            idx2 = convert(Int64,cns_mat[r,2])
            if cns_mat[r,3]==0
                insert!(β,idx1,cns_mat[r,5])
            else
                insert!(β,idx1,cns_mat[r,5]+cns_mat[r,4]*β[idx2])
            end
        end
        # now build the likelihood
        loglike = -sum(-.5 .* ( log(2*π) .+ log(σ^2) .+ ( (y .- X*β)./σ ).^2 ) )
        # more intuitive way? (but JuMP can't use pdf's from Distributions.jl)
        #loglike = -sum( log(1 ./ sqrt(σ^2)) .+ logpdf.(Normal(0,1),(y .- X*β)./sqrt(σ^2)) )
        return loglike
    end

    # run the optimizer for MLE
    # now we need additional information: whether the constraint is "type 1" (set equal to fixed value) or "type 2" (set equal to another parameter)
    # first, set up constraints
    #     Type 1  Restricting one parameter ("parmA") to equal a fixed value
    #     Type 2  Restricting one parameter, parmA, to equal another ("parmB"),
    #             potentially multiplied by some real number q and addd to
    #             some constant m, e.g. parmA = m + q*parmB.
    #
    #   RESTRMAT follows a very specific format. It is an R-by-5 matrix,
    #   where R is the number of restrictions. The role of each of the four
    #   columns is as follows
    #
    #   Column 1  The index of parmA
    #   Column 2  The index of parmB (zero if type 1 restriction)
    #   Column 3  Binary vector w   here 0 indciates a type 1 restriction (parmA
    #               set equal to fixed value) and 1 indicates a type 2
    #               restriction (parmA set equal to parmB)
    #   Column 4  If a type 1 restriction, 0. If     a type 2 restriction, any
    #               real number q such that parmA = q*parmB.
    #   Column 5  If a type 1 restriction, the fixed value. If a type 2
    #                restriction, any real number m such that parmA = m+q*parmB.
    #   NOTE: parmA should always be a later index than parmB
    cns_mat2 = [2 0 0 0 .16;
                4 3 1 2 1]
    svals = vcat(X\y,.5)
    # constraints are treated as data, so take them out of starting values (they get added back in inside the obj function)
    for r=1:size(cns_mat2,1)
        deleteat!(svals, convert(Int64,cns_mat2[r,1]))
    end
    td = TwiceDifferentiable(th -> cns2_reg_mle(th, cns_mat2, X, y), svals; autodiff = :forward)
    θ̂_optim_ad = optimize(td, svals, Newton(linesearch = BackTracking()), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=1))
    θ̂_mle_optim_ad = θ̂_optim_ad.minimizer
    loglikeval = θ̂_optim_ad.minimum
    # evaluate the Hessian at the estimates
    H  = Optim.hessian!(td, θ̂_mle_optim_ad)
    θ̂_mle_optim_ad_se = sqrt.(diag(inv(H)))
    println(θ̂_mle_optim_ad)
    # add back in constraint in both estimates and SEs
    for r in 1:size(cns_mat2,1)
        idx1 = convert(Int64,cns_mat2[r,1])
        idx2 = convert(Int64,cns_mat2[r,2])
        if cns_mat2[r,3]==0
            insert!(θ̂_mle_optim_ad,idx1,cns_mat2[r,5])
            insert!(θ̂_mle_optim_ad_se,idx1,0)
        else
            insert!(θ̂_mle_optim_ad,idx1,cns_mat2[r,5]+cns_mat2[r,4]*θ̂_mle_optim_ad[idx2])
            insert!(θ̂_mle_optim_ad_se,idx1,cns_mat2[r,5]+cns_mat2[r,4]*θ̂_mle_optim_ad_se[idx2]) # this is wrong
        end
        println(θ̂_mle_optim_ad)
    end
    # store results in a data frame
    results.coef_optim_cns2 = vcat(θ̂_mle_optim_ad,-loglikeval)
    results.se_optim_cns2   = vcat(θ̂_mle_optim_ad_se,missing)
    @show results

    # Analytical gradient
    @views @inline function asclogit(bstart::Vector,Y::Array,X::Array,Z::Array,J::Int64,baseAlt::Int64=J,W::Array=ones(length(Y)))

    ## error checking
    @assert ((!isempty(X) || !isempty(Z)) && !isempty(Y))    "You must supply data to the model"
    @assert (ndims(Y)==1 && size(Y,2)==1)                    "Y must be a 1-D Array"
    @assert (minimum(Y)==1 && maximum(Y)==J) "Y should contain integers numbered consecutively from 1 through J"
    if !isempty(X)
        @assert ndims(X)==2          "X must be a 2-dimensional matrix"
        @assert size(X,1)==size(Y,1) "The 1st dimension of X should equal the number of observations in Y"
    end
    if !isempty(Z)
        @assert ndims(Z)==3          "Z must be a 3-dimensional tensor"
        @assert size(Z,1)==size(Y,1) "The 1st dimension of Z should equal the number of observations in Y"
        @assert size(Z,3)==J         "The 3rd dimension of Z should equal the number of choice alternatives"
    end

    K1 = size(X,2)
    K2 = size(Z,2)
    jdx = setdiff(1:J,baseAlt)

    function f(b)
        T = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
        num   = zeros(T,size(Y))
        dem   = zeros(T,size(Y))
        temp  = zeros(T,size(Y))
        numer =  ones(T,size(Y,1),J)
        P     = zeros(T,size(Y,1),J)
        ℓ     =  zero(T)
        b2 = b[K1*(J-1)+1:K1*(J-1)+K2]

        k = 1
        for j in 1:J
            if j != baseAlt
                temp       .= X*b[(k-1)*K1+1:k*K1] .+ (Z[:,:,j].-Z[:,:,baseAlt])*b2
                num        .= (Y.==j).*temp.+num
                dem        .+= exp.(temp)
                numer[:,j] .=  exp.(temp)
                k += 1
            end
        end
        dem.+=1
        P   .=numer./(1 .+ sum(numer;dims=2))

        ℓ = -W'*(num.-log.(dem))
    end

    function g!(G,b)
        T     = promote_type(promote_type(promote_type(eltype(X),eltype(b)),eltype(Z)),eltype(W))
        numer = zeros(T,size(Y,1),J)
        P     = zeros(T,size(Y,1),J)
        numg  = zeros(T,K2)
        demg  = zeros(T,K2)
        b2    = b[K1*(J-1)+1:K1*(J-1)+K2]

        G .= T(0)
        k = 1
        for j in 1:J
            if j != baseAlt
                numer[:,j] .= exp.( X*b[(k-1)*K1+1:k*K1] .+ (Z[:,:,j].-Z[:,:,baseAlt])*b2 )
                k += 1
            end
        end
        P   .=numer./(1 .+ sum(numer;dims=2))

        k = 1
        for j in 1:J
            if j != baseAlt
                G[(k-1)*K1+1:k*K1] .= -X'*(W.*((Y.==j).-P[:,j]))
                k += 1
            end
        end

        for j in 1:J
            if j != baseAlt
                numg .-= (Z[:,:,j].-Z[:,:,baseAlt])'*(W.*(Y.==j))
                demg .-= (Z[:,:,j].-Z[:,:,baseAlt])'*(W.*P[:,j])
            end
        end
        G[K1*(J-1)+1:K1*(J-1)+K2] .= numg.-demg
        return nothing
    end

    td = TwiceDifferentiable(f, g!, bstart, autodiff = :forwarddiff)
    rs = optimize(td, bstart, LBFGS(; linesearch = LineSearches.BackTracking()), Optim.Options(iterations=100_000,g_tol=1e-8,f_tol=1e-8,x_tol=1e-8,show_trace=true))
    β  = Optim.minimizer(rs)
    ℓ  = Optim.minimum(rs)*(-1)
    H  = Optim.hessian!(td, β)
    g  = Optim.gradient!(td, β)
    se = sqrt.(diag(inv(H)))

    return β,se,ℓ,g
    end

    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS4-mixture/nlsw88t.csv"
    dff = CSV.read(HTTP.get(url).body)
    XX = [dff.age dff.white dff.collgrad]
    ZZ = cat(dff.elnwage1, dff.elnwage2, dff.elnwage3, dff.elnwage4,
             dff.elnwage5, dff.elnwage6, dff.elnwage7, dff.elnwage8; dims=3)
    yy = dff.occ_code
    J  = 8
    startvals = [2*rand(7*size(XX,2)).-1; .1]
    β,se,ℓ,g = asclogit(startvals,yy,XX,ZZ,J,J,ones(length(yy)))
    ans = [ .0403744; .2439942; -1.57132; .0433254; .1468556; -2.959103; .1020574; .7473086; -4.12005; .0375628; .6884899; -3.65577; .0204543; -.3584007; -4.376929; .1074636; -.5263738; -6.199197; .1168824; -.2870554; -5.322248; 1.307477]
    dfr = DataFrame(β=β,se=se,answer=ans)
    # time: about 15 seconds
    @show dfr

    @views @inline function mlogit_with_Z(theta, X, Z, y)
        alpha = theta[1:end-1]
        gamma = theta[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end
    ZZ = cat(dff.elnwage1, dff.elnwage2, dff.elnwage3, dff.elnwage4,
             dff.elnwage5, dff.elnwage6, dff.elnwage7, dff.elnwage8; dims=2)
    td = TwiceDifferentiable(theta -> mlogit_with_Z(theta, XX, ZZ, yy), startvals; autodiff = :forward)
    # run the optimizer
    theta_hat_optim_ad = optimize(td, startvals, LBFGS(; linesearch = LineSearches.BackTracking()), Optim.Options(iterations=100_000,g_tol=1e-8,f_tol=1e-8,x_tol=1e-8,show_trace=true))
    theta_hat_mle_ad = theta_hat_optim_ad.minimizer
    # evaluate the Hessian at the estimates
    H  = Optim.hessian!(td, theta_hat_mle_ad)
    theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
    results = DataFrame(coef = theta_hat_mle_ad, se = theta_hat_mle_ad_se)

    # time: about 70 seconds

    return nothing
end

wrapper()


text:
include("../examples/baker.jl")
using LinearAlgebra
using JLD
using SharedArrays
using Distributed
function unstable_sens_gen(s,nSteps)
	u = 2*pi*rand(2)
	u_trj = step(u, s, nSteps)'
	x, y = view(u_trj,:,1), view(u_trj,:,2)
	dJds = 0.
	N = 20
	nSteps = nSteps + 1
	g = zeros(nSteps)
	J = cos.(4*y)
	dXudx1 = zeros(2,nSteps)
	Xu = [sin.(x) zeros(nSteps)]'
	for i = 1:nSteps-1
		x1, x2 = x[i], y[i]
		z = 2.0 + s[1]*cos(x1)
		dzdx = -s[1]*sin(x1)
		dXudx1[:,i] = [cos(x1)/z, 0]
		g[i+1] = g[i]/z - dzdx/z*z
	end
	Xu = Xu[1,:]
	dXudx1 = dXudx1[1,:]
	for n = 1:N
		J_shift = J[n+1:end]
		nJ = length(J_shift)
		g_shift = g[2:nJ+1]
		dXudx1_shift = dXudx1[1:nJ]
		Xu_shift = Xu[1:nJ]
		dJds -= dot(J_shift,dXudx1_shift)/nJ
		dJds -= dot(J_shift,g_shift.*Xu_shift)/nJ
	end
	#@show le, dJds
	return dJds

end

function unstable_sens(s,nSteps)
	u = 2*pi*rand(2)
	u_trj = step(u, s, nSteps)'
	x, y = view(u_trj,:,1), view(u_trj,:,2)
	dJds = 0.
	N = 4
	nSteps = nSteps + 1
	g = zeros(nSteps)
	J = cos.(4*y)
	dXudx1 = zeros(2,nSteps)
	Xu = [sin.(x) zeros(nSteps)]'
	for i = 1:nSteps-1
		x1, x2 = x[i], y[i]
		z = 2.0 + s[1]*cos(x1)
		dzdx = -s[1]*sin(x1)
		dXudx1[:,i] = [cos(x1)/z, 0]
		g[i+1] = g[i]/z - dzdx/z*z
	end
	Xu = Xu[1,:]
	dXudx1 = dXudx1[1,:]
	for n = 1:N
		J_shift = J[n+1:end]
		nJ = length(J_shift)
		g_shift = g[2:nJ+1]
		dXudx1_shift = dXudx1[1:nJ]
		Xu_shift = Xu[1:nJ]
		dJds -= dot(J_shift,dXudx1_shift)/nJ
		dJds -= dot(J_shift,g_shift.*Xu_shift)/nJ
	end
	#@show le, dJds
	return dJds

end
function get_unstable_sens(s1)
	p = 4
	nSteps = 500000
	s = zeros(p)
	# J = cos(4y)
	n_exps = size(s1)[1]
	dJds = zeros(n_exps)
	n_rep = 160
	dJds_proc = SharedArray{Float64}(n_rep)
	for k=1:n_exps
		s[1] = s1[k]
		dJds_proc .= 0.
		t = @distributed for j=1:n_rep
			dJds_proc[j] = unstable_sens(s,nSteps)/n_rep
		end
		wait(t)
		dJds[k] = sum(dJds_proc)
		@show dJds[k]
	end
	save("../data/unstable_sens/dJds1_test.jld", "s1",
	     s1, "dJds", dJds)
end

ext: jl

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
"""
    simplex_projection_and_support(z)

Compute the Euclidean projection onto the probability simplex and the set of indices where it is nonzero.

See <https://arxiv.org/abs/1602.02068> for details.
"""
function simplex_projection_and_support(z::AbstractVector{R}) where {R<:Real}
    d = length(z)
    z_sorted = sort(z; rev=true)
    z_sorted_cumsum = cumsum(z_sorted)
    k = maximum(j for j in 1:d if (1 + j * z_sorted[j]) > z_sorted_cumsum[j])
    τ = (z_sorted_cumsum[k] - 1) / k
    p = Vector{R}(undef, d)
    s = Vector{Int}(undef, d)
    for i in 1:d
        p[i] = max(z[i] - τ, zero(R))
        s[i] = Int(!iszero(p[i]))
    end
    return p, s
end;

"""
    simplex_projection(z)

Compute the Euclidean projection onto the probability simplex.
"""
function simplex_projection(z::AbstractVector{<:Real})
    p, _ = simplex_projection_and_support(z)
    return p
end;

"""
    rrule(::typeof(simplex_projection), z)

Custom reverse rule for [`simplex_projection`](@ref) which bypasses the sorting step.

See <https://arxiv.org/abs/1602.02068> for details.
"""
function ChainRulesCore.rrule(::typeof(simplex_projection), z::AbstractVector{<:Real})
    p, s = simplex_projection_and_support(z)
    S = sum(s)
    function simplex_projection_pullback(dp)
        vjp = s .* (dp .- (dp's) / S)
        return (NoTangent(), vjp)
    end
    return p, simplex_projection_pullback
end;


text:
module EOS

using AbInitioSoftwareBase: load
using Comonicon: @cast
using EquationOfStateRecipes
using EquationsOfStateOfSolids: Parameters, EquationOfStateOfSolids, EnergyEquation
using Express.EquationOfStateWorkflow: Scf, VcOptim, FitEos, buildjob
using Express.EquationOfStateWorkflow.Recipes: buildworkflow
using QuantumESPRESSOExpress
using Serialization: deserialize
using SimpleWorkflows: run!, getresult
import Plots
# using ..ExpressCommands: @load_plugin

"""
Fit an equation of state from `cfgfile` for calculation `calc`.

# Arguments
- `calc`: the calculation type. Acceptable options are "scf" for self-consistent calculation and "optim" for structure optimizations.
- `cfgfile`: the file to be printed. Available extensions are `.jls`, `.json`, `.yaml`, `.yml` or `.toml`.
"""
@cast function fit(calc, cfgfile)
    calc = lowercase(calc)
    if calc == "scf"
        T = Scf
    elseif calc == "optim"
        T = VcOptim
    else
        throw(ArgumentError("unrecognized calculation type `$calc`!"))
    end
    job = buildjob(FitEos{T}(), cfgfile)
    run!(job)
    wait(job)
    display(something(getresult(job)))
end

"""
Plot an equation of state from `file` to file "out".

# Arguments
- `file`: a `.jls` file that saves the equation of state.
"""
@cast function plot(file)
    data = deserialize(file)
    for (key, value) in data
        if value isa Parameters
            eos = EnergyEquation(value)
        elseif data isa EquationOfStateOfSolids
            eos = value
        else
            error("unsupported type $(typeof(value))!")
        end
        Plots.plot(eos)
        Plots.savefig("$key.pdf")
    end
end

"""
Run a `config` file, better with absolute path.

# Arguments
- `config`: the file to be run. Available extensions are `.json`, `.yaml`, `.yml` or `.toml`.
"""
@cast function run(config)
    # @load_plugin
    workflow = buildworkflow(config)
    dict = load(config)
    run!(workflow; saveas = dict["save"]["status"])
    return workflow
end

end


text:
const version = v"0.1.0"

# Auth
const TRUST_ON_FIRST_USE = 0  # Deprecated
const TRUST_SIGNED_CERTIFICATES = 1  # Deprecated
const TRUST_ALL_CERTIFICATES = 2
const TRUST_CUSTOM_CA_SIGNED_CERTIFICATES = 3
const TRUST_SYSTEM_CA_SIGNED_CERTIFICATES = 4
const TRUST_DEFAULT = TRUST_ALL_CERTIFICATES

# Connection Pool Management
const INFINITE = -1
const DEFAULT_MAX_CONNECTION_LIFETIME = 3600  # 1h
const DEFAULT_MAX_CONNECTION_POOL_SIZE = 100
const DEFAULT_CONNECTION_TIMEOUT = 5.0  # 5s

# Connection Settings
const DEFAULT_CONNECTION_ACQUISITION_TIMEOUT = 60  # 1m

# Routing settings
const DEFAULT_MAX_RETRY_TIME = 30.0  # 30s

const LOAD_BALANCING_STRATEGY_LEAST_CONNECTED = 0
const LOAD_BALANCING_STRATEGY_ROUND_ROBIN = 1
const DEFAULT_LOAD_BALANCING_STRATEGY = LOAD_BALANCING_STRATEGY_LEAST_CONNECTED

# Client name
const DEFAULT_USER_AGENT = "juliabolt/$(string(version)) Julia/$(string(VERSION)) ($(string(Sys.KERNEL)))"

const default_config = Dict(
    "auth" => nothing,  # provide your own authentication token such as {"username", "password"}
    "encrypted" => nothing,  # default to have encryption enabled if ssl is available on your platform
    "trust" => TRUST_DEFAULT,
    "der_encoded_server_certificate" => nothing,

    "user_agent" => DEFAULT_USER_AGENT,

    # Connection pool management
    "max_connection_lifetime" => DEFAULT_MAX_CONNECTION_LIFETIME,
    "max_connection_pool_size" => DEFAULT_MAX_CONNECTION_POOL_SIZE,
    "connection_acquisition_timeout" => DEFAULT_CONNECTION_ACQUISITION_TIMEOUT,

    # Connection settings:
    "connection_timeout" => DEFAULT_CONNECTION_TIMEOUT,
    "keep_alive" => nothing,

    # Routing settings:
    "max_retry_time" => DEFAULT_MAX_RETRY_TIME,
    "load_balancing_strategy" => DEFAULT_LOAD_BALANCING_STRATEGY,
)


text:
using ..CryptoGroups: static

struct FP{P} <: PrimeField
    x::BigInt # Instead of BigInt I could use BitIntegers here. Need to check performance of powermod...
    FP{P}(x::Integer) where P = new{P}(BigInt(x))
end

specialize(::Type{FP}, p::Integer) = FP{static(p)}

value(x::FP) = x.x
modulus(::FP{P}) where P = BigInt(P)
modulus(::Type{FP{P}}) where P = BigInt(P)

# Why order and modulus is the same thing here?
order(::FP{P}) where P = BigInt(P)
order(::Type{FP{P}}) where P = BigInt(P)

function trimnumber(x::String)
    if length(x) < 30
        return x
    else
        return x[1:10] * "..." * x[end-10:end]
    end
end

trimnumber(x::Integer)= trimnumber(string(x))

modulus_str(x::BigInt) = trimnumber(x)
modulus_str(::Type{F}) where F <: FP = modulus_str(modulus(F))

Base.show(io::IO, ::Type{FP}) = print(io, "FP")
Base.show(io::IO, ::Type{F}) where F <: FP = print(io, "FP/$(modulus_str(F))")

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

### StaticBitVEctor would be a better fit here

struct Reducer{N}
    len::Int
    chunks::NTuple{N, UInt64}
end

Reducer(x::BitVector) = Reducer(x.len, Tuple(x.chunks))

Reducer(x::Vector{Int}) = Reducer(BitVector(i in x for i in 0:maximum(x)))

Base.length(x::Reducer) = x.len

bitlength(R::Reducer) = length(R) - 1

function Base.convert(::Type{BitVector}, x::Reducer)

    a = BitVector(undef, x.len)
    a.chunks = UInt64[x.chunks...]

    return a
end

function xstr(i::Int)
    if i == 0
        return "1"
    elseif i == 1
        return "X"
    else
        return "X^$i"
    end
end

get_powers(x::BitVector) = (0:x.len - 1)[x]
get_powers(x::Reducer) = get_powers(convert(BitVector, x))

function print_poly(io::IO, x::BitVector)
    n = reverse(get_powers(x)) # Best to visualize dominant first
    write(io, "$(xstr(n[1]))")
    for i in n[2:end]
        write(io, " + $(xstr(i))")
    end
end

print_poly(io::IO, x::Reducer) = print_poly(io, convert(BitVector, x))
print_poly(x::Any) = print_poly(Base.stdout, x)

Base.show(io::IO, x::Reducer) = print_poly(io, x)

@enum Endian big little

struct F2PB{R} <: BinaryField ### R is reducer
    x::BitVector

    function F2PB{R}(x::BitVector, e::Endian) where R

        @assert length(x) == length(R) - 1

        if e == big
            new(x)
        elseif e == little
            new(reverse(x))
        end
    end

    F2PB{R}(x::BitVector; e::Endian = big) where R = F2PB{R}(x, e)

    F2PB(f) = F2PB{Reducer(f)} # Perhaps I could get rid of Reducer in favour of a StaticBitVector
    F2PB(f, x) = F2PB(f)(x)
end

F2PB(x::F2PB) = x

bitlength(::Type{F2PB{R}}) where R = length(R) - 1
bitlength(::F) where F <: F2PB = bitlength(F)

Base.convert(::Type{F}, x::BitVector) where F <: F2PB = F(x, little)

# This function is necessary in differetn contexts so I could keep it as is.
tobits(x::F2PB) = reverse(x.x)

Base.zero(::Type{F2PB{R}}) where R = F2PB{R}(BitVector(0 for i in 1:length(R)-1))
Base.one(::Type{F2PB{R}}) where R = F2PB{R}(BitVector((1, (0 for i in 2:length(R)-1 )...)))

Base.:+(x::F, y::F) where F <: F2PB = F(xor.(x.x, y.x))

order(::Type{F}) where F <: F2PB = BigInt(2)^bitlength(F) - 1
order(x::F) where F <: F2PB = order(F)

_order(x::BitVector) = findlast(x->x==1, x) - 1
_rem(x::BitVector) = x[1:findlast(x->x==1, x)-1]

function red!(c::BitVector, f::BitVector)

    m = _order(f)

    r = _rem(f)
    append!(r, BitVector(0 for i in 1:m))

    s = 2m - length(c)

    append!(c, BitVector(0 for i in 1:s))

    for i in length(c):-1:m+1 ### iteration may be reverse order
        if c[i] == true
            ri = r >> ((i - 1) - m )
            c = xor.(c, ri)
        end
    end

    return c[1:m]
end

function mul(a::BitVector, b::BitVector)

    a = copy(a)

    n = length(a)

    append!(a, BitVector(0 for i in 1:n))

    s = zero(a)

    for i in 1:n
        if b[i]
            s .⊻= a >> (i - 1)
        end
    end

    return s
end

mul(a::BitVector, b::BitVector, f::BitVector) = red!(mul(a, b), f)
Base.:*(x::F2PB{R}, y::F2PB{R}) where R = F2PB{R}(mul(x.x, y.x, convert(BitVector, R)))
Base.:(==)(x::F, y::F) where F <: F2PB = x.x == y.x

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

struct F2GNB{N, T} <: BinaryField

    x::BitVector

    function F2GNB{N, T}(x::BitVector) where {N, T}
        @assert length(x) == N
        new(x)
    end

end

Base.convert(::Type{F2GNB{N, T}}, bits::BitVector) where {N, T} = F2GNB{N, T}(bits)

tobits(x::F2GNB) = x.x

bitlength(::Type{F2GNB{N, T}}) where {N, T} = N
bitlength(::F) where F <: F2GNB = bitlength(F)

order(::Type{F}) where F <: F2GNB = BigInt(2)^bitlength(F) - 1
order(x::F) where F <: F2GNB = order(F)

Base.zero(::Type{F2GNB{N, T}}) where {N, T} = F2GNB{N, T}(BitVector(false for i in 1:N))
Base.one(::Type{F2GNB{N, T}}) where {N, T} = F2GNB{N, T}(BitVector(true for i in 1:N))

square(x::F) where F <: F2GNB = F(circshift(x.x, 1))

Base.:+(x::F, y::F) where F <: F2GNB = F(xor.(x.x, y.x))
Base.:(==)(x::F, y::F) where F <: F2GNB = x.x == y.x

# Improving performace of this function is top priority
# Random access order with Fvec is likelly making a lot of chache misses

function first_coordinate_product(a::BitVector, b::BitVector, T::P, F::Vector{P}) where P <: Integer

    m = length(a)

    p = T*m + 1

    J = 0

    if mod(T, 2) != 0

        # In the spec it is given that summation index goes to m/2
        # If m/2 would be performed as ceiling
        # The following code would ask for an element like b[m + 1]
        # Alternativelly a leading bit could be added then to make araray even
        for k in 1:div(m, 2)
            J ⊻= a[k] * b[div(m, 2) + k]
            J ⊻= a[div(m, 2) + k] * b[k]
        end
    end

    c0 = J

    for k in 1:(p - 2)
        c0 ⊻= a[F[k + 1] + 1] * b[F[p - k] + 1]
    end

    return c0

end

# Is a dublicate, I could import it from Specs when dust settles
function compute_integer_order(g::T, p::T) where T <: Integer

    @assert 1 < g < p

    b = g

    j = 1 ### Could there be an error?

    while b > 1
        b = mod(g*b, p)
        j += 1
    end

    return j
end

function construct_integer_order_prime(T, p)

    for g in 2:(p-1)

        k = compute_integer_order(g, p)

        if mod(k, T) == 0
            u = powermod(g, div(k, T), p)
            return u
        end
    end

    error("No success")
end

function compute_F(m::Int, T::Int)

    p = T*m + 1

    u = construct_integer_order_prime(T, p)

    F = Vector{Int}(undef, p)

    w = 1

    for j in 0:(T - 1)

        n = w

        for i in 0:(m - 1)
            F[n] = i
            n = mod(2n, p)
        end

        w = mod(u*w, p)
    end

    return F
end

function mul_gnb(a::BitVector, b::BitVector, T::Int)

    u = copy(a)
    v = copy(b)

    m = length(a)

    F = compute_F(m, T)

    c = zero(u)

    # botleneck here
    for k in 0:(m - 1)
        c[k+1] = first_coordinate_product(u, v, T, F)
        circshift!(u, -1)
        circshift!(v, -1)
    end

    return c
end

Base.:*(x::F2GNB{N, T}, y::F2GNB{N, T}) where {N, T} = F2GNB{N, T}(mul_gnb(x.x, y.x, T))


text:
function LLop(basis,quad)

    dim = dimension(basis)
    sdim = symmetric_tensor_dimension(dim)

    return mass_matrix(basis,quad,-1.0,sdim)
end

function LLop(basis,quad,cellmap)
    return determinant_jacobian(cellmap)*LLop(basis,quad)
end

function LUop(basis,grad_basis,quad,ED,cellmap,NF)

    @assert length(ED) > 0

    dim,sdim = size(ED[1])
    @assert all([size(E) == (dim,sdim) for E in ED])

    ALU = zeros(sdim*NF,dim*NF)
    invjac = inverse_jacobian(cellmap)
    detjac = determinant_jacobian(cellmap)

    for (p,w) in quad
        vals = basis(p)
        grads = grad_basis(p)
        for k = 1:length(ED)
            gradsk = invjac[k]*grads[:,k]
            Mk = make_row_matrix(gradsk,ED[k])
            N = interpolation_matrix(vals,dim)

            ALU .+= Mk'*N*detjac*w
        end
    end

    return ALU
end

function LUop(basis,quad,Dhalf,cellmap)
    dim = dimension(basis)
    NF = number_of_basis_functions(basis)
    Ek = vec_to_symm_mat_converter(dim)
    ED = [E'*Dhalf for E in Ek]
    return LUop(basis,x->gradient(basis,x),quad,ED,cellmap,NF)
end

function UUop(basis,facequads,facemaps,stabilization,cellmap)
    dim = dimension(basis)
    scale = face_determinant_jacobian(cellmap)
    return stabilization*mass_matrix_on_boundary(basis,facequads,facemaps,
        scale,dim)
end

function UUop(basis,facequads,facemaps,iquad,inormals,imap,stabilization,
    cellmap,ndofs,NF)

    facescale = face_determinant_jacobian(cellmap)
    iscale = scale_area(cellmap,inormals)
    return stabilization*mass_matrix_on_boundary(basis,facequads,facemaps,
        facescale,iquad,imap,iscale,ndofs,NF)

end

function UUop(basis,facequads,facemaps,iquad,inormals,imap,
    stabilization,cellmap)

    dim = dimension(basis)
    nf = number_of_basis_functions(basis)

    return UUop(basis,facequads,facemaps,iquad,inormals,imap,stabilization,
        cellmap,dim,nf)
end

function local_operator(LL,LU,UU)
    lop = [LL   LU
           LU'  UU]
    return lop
end

function local_operator(basis,vquad,facequad,facemaps,Dhalf,
    stabilization,cellmap)

    LL = LLop(basis,vquad,cellmap)
    LU = LUop(basis,vquad,Dhalf,cellmap)
    UU = UUop(basis,facequad,facemaps,stabilization,cellmap)
    return local_operator(LL,LU,UU)
end

function local_operator(basis,vquad,facequads,facemaps,iquad,imap,inormals,
    Dhalf,stabilization,cellmap)

    LL = LLop(basis,vquad,cellmap)
    LU = LUop(basis,vquad,Dhalf,cellmap)
    UU = UUop(basis,facequads,facemaps,iquad,inormals,imap,
        stabilization,cellmap)
    return local_operator(LL,LU,UU)
end

function body_force_linear_form(vbasis,vquad,bodyforce,cellmap)
    dim = dimension(vbasis)
    NF = number_of_basis_functions(vbasis)
    sdim = symmetric_tensor_dimension(dim)
    @assert length(bodyforce) == dim
    rhs = determinant_jacobian(cellmap)*linear_form(bodyforce,vbasis,vquad)
    return vcat(zeros(sdim*NF),rhs)
end


text:
"""
```
chull = convexhull(img)
```
Computes the convex hull of a binary image and returns the vertices of convex hull as a CartesianIndex array.
"""
function convexhull(img::AbstractArray{T, 2}) where T<:Union{Bool,Gray{Bool}}

    function getboundarypoints(img)
        points = CartesianIndex{2}[]
        for j = axes(img, 2)
            v = Base.view(img, :, j)
            i1 = findfirst(v)
            if i1 != nothing
                i2 = findlast(v)
                push!(points, CartesianIndex(i1, j))
                if i1 != i2
                    push!(points, CartesianIndex(i2, j))
                end
            end
        end
        points
    end

    function right_oriented(ref, a, b)
        (a[2]-ref[2])*(b[1]-ref[1])-(b[2]-ref[2])*(a[1]-ref[1])<0
    end

    function collinear(a, b, c)
        (a[2]-c[2])*(b[1]-c[1])-(b[2]-c[2])*(a[1]-c[1])==0
    end

    dist2(a, b) = sum(abs2, (a-b).I)

    function angularsort!(points, ref, img)
        last_point = ref

        for i in eachindex(points)
            if points[i]==last_point
                deleteat!(points, i)
                break
            end
        end

        sort!(points, lt=(a, b) -> right_oriented(last_point, a, b))

        i=0
        while i<=length(points)
            i=i+1
            if i+1>length(points)
                break
            end

            if collinear(last_point, points[i], points[i+1])
                if dist2(last_point, points[i]) < dist2(last_point, points[i+1])
                    deleteat!(points, i)
                else
                    deleteat!(points, i+1)
                end
                i=i-1
            end
        end
    end

    # Used Graham scan algorithm

    points=getboundarypoints(img)
    last_point=CartesianIndex(map(r->first(r)-1, axes(img)))
    for point in points
        if point[2]>last_point[2] || (point[2]==last_point[2] && point[1]>last_point[1])
            last_point=point
        end
    end
    angularsort!(points, last_point, img)

    n_points=length(points)

    if n_points<3
        error("Not enough points to compute convex hull.")
    end

    convex_hull=CartesianIndex{2}[]
    push!(convex_hull, last_point)
    push!(convex_hull, points[1])
    push!(convex_hull, points[2])

    for i in 3:n_points
        while (right_oriented(convex_hull[end],convex_hull[end-1], points[i]) || collinear(convex_hull[end],convex_hull[end-1], points[i]))
            pop!(convex_hull)
        end
        push!(convex_hull, points[i])
    end

    return convex_hull
end


text:
"""
    MultivariateGammaConvolution(α,θ)

Constructs a distribution that corresponds to the convolutions of MutlivariateGamma(α[i],θ[i,:]) distributions.
The distribution can then be used through several methods, following the Distributions.jl standard

The pdf and cdf and not yet coded.
The code is type stable and handles any <:Real types, given by the parameters.

Fitting th edistribution is not yet possible.

# Examples
```julia-repl
julia> dist = MultivariateGammaConvolution([2,3,1],[3 0;4 0;0 1]);
julia> sample = zeros(Float64,(2,10));
julia> Random.rand!(dist,sample);

```
"""
struct MultivariateGammaConvolution{T<:Real, V<:AbstractVector{T}, M<: AbstractMatrix{T}} <: Distributions.ContinuousMultivariateDistribution where T
    α::V
    θ::M
    constants::Int
end

n(x::MultivariateGammaConvolution) = length(x.α)

function Base.show(io::IO, m::MultivariateGammaConvolution)
    println("Multivariate Gamma Convolutions with parametrisation:")
    display([m.α m.θ])
end

#### Constructor :
function MultivariateGammaConvolution(α::AbstractVector{T1},θ::AbstractMatrix{T2}) where {T1 <: Real, T2<:Real}

     @assert(length(α) == size(θ,1),"You must provide α with same length as the number of rows in θ")
     # Propote eltypes :
     T = Base.promote_eltype(α,θ,[1.0])
     α = T.(α)
     θ = T.(θ) # the size of theta is (n,d)

     # Remove alphas that are non-positives or thetas that are all zero.
     if any(α .< eps(T))
         are_pos = α .> eps(T) | [all(θ[i,:] .> 0) for i in 1:size(θ,1)]
         α = α[are_pos]
         θ = θ[are_pos,:]
     end

     if length(α) == 1
         return MultivariateGamma(α,θ[:,1])
     end
     # sorting by alphas :
     θ =θ[sortperm(α),:]
     sort!(α)
     return MultivariateGammaConvolution(α,θ, 1)
end
MultivariateGammaConvolution(α::T1,θ::AbstractVector{T2}) where {T1<:Real, T2<:Real} = MultivariateGamma(promote(α,θ)...)
MultivariateGammaConvolution(α::AbstractVector{T1},θ::AbstractVector{T2}) where {T1<:Real, T2<:Real} = UnivariateGammaConvolution(promote(α,θ)...)

#### eltype, length, support
Base.eltype(d::MultivariateGammaConvolution) = typeof(d.α)
Base.length(d::MultivariateGammaConvolution) = size(d.θ,2)
function Distributions.insupport(d::MultivariateGammaConvolution,x::AbstractVector{T}) where {T <: Real}
    return all(x > T(0))
end

#### Sampling
struct MGCSPL <: Distributions.Sampleable{Distributions.Multivariate,Distributions.Continuous}
    Γs::AbstractVector{MGSPL}
end
Distributions.sampler(d::MultivariateGammaConvolution) = MGCSPL([Distributions.sampler(MultivariateGamma(d.α[i],d.θ[i,:])) for i in 1:length(d.α)])
Base.length(s::MGCSPL) = length(s.Γs[1].θ)
function Distributions._rand!(rng::Distributions.AbstractRNG, s::MGCSPL, x::AbstractVector{T}) where T<:Real
    x .= T(0)
    for i in 1:length(s.Γs)
        x .+= rand(rng,s.Γs[i])
    end
end

#### Logpdf
function Distributions._logpdf(d::MultivariateGammaConvolution{T}, x::AbstractArray) where T<:Real
    if !insupport(d,x)
        return T(-Inf)
    end
    error("Not implemented yet.")
end

# mean(d::MultivariateGamma) = d.α .* d.θ
# var(d::MultivariateGamma) = d.α .* d.θ .* d.θ
# Statistics.cov(d::MultivariateGamma) = d.α .* d.θ[na,:] .* d.θ[:,na]

# Okay this should be enough to obtain a pdf

# dist = MultivariateGammaConvolution([2,3,1],[3 0;4 0;0 1])
# sample = zeros(Float64,(2,10))
# import Random
# Random.rand!(dist,sample)
# display(sample)
#
#
# display(Distributions.pdf(dist,sample))
# display(Distributions.pdf(Distributions.Gamma(dist.α,dist.θ[1]),sample[1,:]))


text:
struct DataBinding

    id    :: Symbol

    here  :: DataUnit
    there :: DataUnit

    # views
    here_pull_view   :: Any
    there_pull_view  :: Any

    here_push_view   :: Any
    there_push_view  :: Any

    # pull from there
    # push from here

    function DataBinding(

        id          :: Symbol,

        here        :: DataUnit,
        there       :: DataUnit,

        here_pull_yrng        :: Any,
        there_pull_yrng       :: Any,

        here_push_yrng        :: Any,
        there_push_yrng       :: Any,
    )

        here_pull_view  = view( here.data, :, :,  here_pull_yrng)
        there_pull_view = view(there.data, :, :, there_pull_yrng)

        here_push_view  = view( here.data, :, :,  here_push_yrng)
        here_push_view  = view(there.data, :, :, there_push_yrng)

        if length(here_pull_view) != length(there_pull_view)
            println(size(here_pull_view), "; " , size(there_pull_view))
            throw(ErrorException("Range mismatch in pull view. Key: " * string(id)))
        end

        if length(here_push_view) != length(there_push_view)
            println(size(here_push_view), "; " , size(there_push_view))
            throw(ErrorException("Range mismatch in push view. Key: " * string(id)))
        end

        return new(
            id,
            here,
            there,
            here_pull_view,
            there_pull_view,
            here_push_view,
            there_push_view,
        )

    end

end

struct DataExchanger

    groups   :: Dict
    bindings :: Dict

    function DataExchanger(
        group_labels :: AbstractArray{Symbol} = Array{Symbol}(undef, 0),
    )

        if ! (:ALL in group_labels)
            push!(group_labels, :ALL)
        end

        groups = Dict()
        for label in group_labels
            groups[label] = Array{Symbol}(undef, 0)
        end

        bindings = Dict()

        return new(
            groups,
            bindings,
        )
    end

end

function createBinding!(
    data_exchanger   :: DataExchanger,
    id               :: Symbol,
    here             :: DataUnit,
    there            :: DataUnit,
    here_pull_yrng        :: Any,
    there_pull_yrng       :: Any,
    here_push_yrng        :: Any,
    there_push_yrng       :: Any;
)

    db = DataBinding(
        id,
        here,
        there,
        here_pull_yrng,
        there_pull_yrng,
        here_push_yrng,
        there_push_yrng,
    )

    if ! haskey(data_exchanger.bindings, id)
        data_exchanger.bindings[id] = db
    else
        throw(ErrorException("Binding id `" * string(id) * "` already existed."))
    end

end

function addToGroup!(
    data_exchanger :: DataExchanger,
    id    :: Symbol,
    label :: Symbol,
)

    if ! haskey(data_exchanger.groups, label)
        println("Group `" * string(label) * "` does not exist. Create one.")
        data_exchanger.groups[label] = Array{Symbol}(undef, 0)
    end

    if id in data_exchanger.groups[label]
        throw(ErrorException("Binding id `" * string(id) * "` already in group `" * string(label) * "`."))
    end

    push!(data_exchanger.groups[label], id)

end

function syncData!(
    data_exchanger :: DataExchanger,
    group_label    :: Symbol,
    direction      :: Symbol,
)

    dbs = data_exchanger.groups[group_label]

    if direction == :PUSH
        for db_label in dbs
            db = data_exchanger.bindings[db_label]

            #=
            if db_label == :Φ
                db.here_push_view .= rand(size(db.here_push_view)...)
                println("Before SYNC: ")
                println("SUM here : ",  SUM(db.here_push_view))
                println("SUM there: ", SUM(db.there_push_view))
            end
            =#

            db.there_push_view .= db.here_push_view

            #=
            if db_label == :Φ

                for i=1:144, j=1:90
                    db.there_push_view[i, j] = db.here_push_view[i, j]
                end

                println("db: id: ", db.here.id, "; ", db.there.id, "; size: ", size(db.there_push_view), "; sizehere: ", size(db.here_push_view))
                println("SUM here : ",  SUM(db.here_push_view))
                println("SUM there: ", SUM(db.there_push_view))
            end
            =#
        end
    elseif direction == :PULL
        for db_label in dbs

            db = data_exchanger.bindings[db_label]
            db.here_pull_view .= db.there_pull_view

#            if db.here.id == :SWFLX
#                println("db: id: ", db.here.id, "; ", db.there.id, "; size: ", size(db.there_push_view), "; sizehere: ", size(db.here_push_view))
#                println("sum: ", sum(db.there.odata))
#            end
        end
    else
        throw(ErrorException("Unknown direction keyword: " * string(direction)))
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# ---
# title: 675. Cut Off Trees for Golf Event
# id: problem675
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/cut-off-trees-for-golf-event/description/>
# hidden: true
# ---
#
# You are asked to cut off all the trees in a forest for a golf event. The
# forest is represented as an `m x n` matrix. In this matrix:
#
#   * `0` means the cell cannot be walked through.
#   * `1` represents an empty cell that can be walked through.
#   * A number greater than `1` represents a tree in a cell that can be walked through, and this number is the tree's height.
#
# In one step, you can walk in any of the four directions: north, east, south,
# and west. If you are standing in a cell with a tree, you can choose whether to
# cut it off.
#
# You must cut off the trees in order from shortest to tallest. When you cut off
# a tree, the value at its cell becomes `1` (an empty cell).
#
# Starting from the point `(0, 0)`, return _the minimum steps you need to walk
# to cut off all the trees_. If you cannot cut off all the trees, return `-1`.
#
# You are guaranteed that no two trees have the same height, and there is at
# least one tree needs to be cut off.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/26/trees1.jpg)
#
#
#
#     Input: forest = [[1,2,3],[0,0,4],[7,6,5]]
#     Output: 6
#     Explanation: Following the path above allows you to cut off the trees from shortest to tallest in 6 steps.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/26/trees2.jpg)
#
#
#
#     Input: forest = [[1,2,3],[0,0,0],[7,6,5]]
#     Output: -1
#     Explanation: The trees in the bottom row cannot be accessed as the middle row is blocked.
#
#
# **Example 3:**
#
#
#
#     Input: forest = [[2,3,4],[0,0,5],[8,7,6]]
#     Output: 6
#     Explanation: You can follow the same path as Example 1 to cut off all the trees.
#     Note that you can cut off the first tree at (0, 0) before making any steps.
#
#
#
#
# **Constraints:**
#
#   * `m == forest.length`
#   * `n == forest[i].length`
#   * `1 <= m, n <= 50`
#   * `0 <= forest[i][j] <= 109`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using impedance
using Test

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

max_stars_repo_name: nirmitdeshpande/EISModeling

max_issues_repo_name: nirmitdeshpande/EISModeling

max_forks_repo_name: nirmitdeshpande/EISModeling


text:
# Exercise 1
#various rules for differentiation

struct DualNumber{T} <: Real
    val::T
    ϵ::T
end
import Base: +, *, -, ^, /, exp
+(x::DualNumber, y::DualNumber) = DualNumber(x.val + y.val, x.ϵ + y.ϵ)  # dual addition
+(x::DualNumber, a::Number) = DualNumber(x.val + a, x.ϵ)  # i.e. scalar addition, not dual
+(a::Number, x::DualNumber) = DualNumber(x.val + a, x.ϵ)  # i.e. scalar addition, not dual
-(x::DualNumber, y::DualNumber) = DualNumber(x.val - y.val, x.ϵ - y.ϵ)  # dual subtraction
-(x::DualNumber, a::Number) = DualNumber(x.val - a, x.ϵ)  # i.e. scalar subtraction, not dual
-(a::Number, x::DualNumber) = DualNumber(x.val - a, x.ϵ)  # i.e. scalar subtraction, not dual
/(x::DualNumber, y::DualNumber) = DualNumber(x.val/y.val, (y.val*x.ϵ - x.val*y.ϵ)/y.val^2) #dual quotient rule
*(x::DualNumber, y::DualNumber) = DualNumber(x.val*y.val, x.val*y.ϵ + y.val*x.ϵ ) #dual product rule
*(x::DualNumber, a::Number) = DualNumber(a*x.val, a*x.ϵ ) #scalar multiplicaiton
^(a::Number, x::DualNumber) = DualNumber(x.val^a, a*x.val^(a-1)*x.ϵ)

f(x) = x^2.0

x = DualNumber(1.0, 1.0)  # x -> 2.0 + 1.0\epsilon
y = DualNumber(3.0, 1.0)  # i.e. y = 3.0, no derivative

# seeded calculates both the function and the d/dx gradient!
f(x)

max_stars_repo_name: shanemcmiken/quantecon-notebooks-julia

max_issues_repo_name: shanemcmiken/quantecon-notebooks-julia

max_forks_repo_name: shanemcmiken/quantecon-notebooks-julia


text:
using Revise
using SeniorProjectAssignments
using DataFrames
using CSV

using SeniorProjectAssignments: match

using Random
# Random.seed(146)

survey_dataframe = CSV.read("survey_012.csv", DataFrame, header=1, datarow=4)
project_dataframe = CSV.read("projects_012.csv", DataFrame)

projects = process_projects(project_dataframe, min=10, max=12)
@show pnames = [p.id for p in projects]

# Random.seed!(146)
# function map_roles(roles)
#     map(roles) do rating
#         round(rating/10.0)
#     end
# end

Random.seed!(147)
map_roles(roles) = map(rating->1.0, roles)
#
students, groups = process_survey(survey_dataframe, pnames, map_roles=map_roles)

# Random.seed!(147)
# students, groups = process_survey(survey_dataframe, pnames)

CSV.write("processed_students.csv", convert(DataFrame, students))

@show groups

sdf = convert(DataFrame, students)

force = ["Andy Benham"=>"Prof. Nerem: Drone Bathymetry"]

m = match(students, projects, groups, force=force)
sort!(m, :project)
CSV.write("output_012.csv", m)
display(m)
display([r=>sum(skipmissing(m[:,:rank]).==r) for r in sort(unique(m[:,:rank]))])


text:
struct ModalInstance <: AbstractVector{SubDataFrame}
    # worlds::Set{<:AbstractOntology} # TODO add set of worlds, for the semantics (RBC)
    rows::Vector{DataFrameRow}
end

ModalInstance() = ModalInstance(DataFrameRow[])
DataFrames.nrow(::Array{DataFrameRow,1}) = 1

function Base.show(io::IO, ::MIME"text/plain", mi::ModalInstance)
    println(io, "Modal instance with $(length(mi.rows)) frames")
    for i in 1:length(mi.rows)
        println(io, "Frame #$i")
        println(io, mi.rows[i])
    end
end

Base.size(mi::ModalInstance) = (nrow(mi.rows),)
Base.getindex(mi::ModalInstance, i::Int) = mi.rows[i]

mutable struct ModalFrame <: AbstractVector{SubDataFrame}
    # name::String # TODO add name?
    dimension::Int  # (0=static, 1=timeseries, 2=images, ecc.)
    data::DataFrame
    views::Vector{DataFrameRow} # intraframe instances (i.e., first, second, .., series)
end

dimension(mf::ModalFrame) = mf.dimension

@inline function ModalFrame(data::DataFrame)
    colnames = names(data)

    @assert length(colnames) > 0 && nrow(data) > 0 "There must be at least one column and one instance in data to compute its dimension."

    dimension = ndims(data[1, Symbol(colnames[1])])

    for colname in colnames
        curr_dimension = ndims(data[1, Symbol(colname)])
        @assert dimension == curr_dimension "Each attribute must have the same dimension."
    end

    # if dimension > 0
    #     typ = eltype(data[1, Symbol(colnams[1])])
    #     for colname in colnames
    #         curr_typ = eltype(data[1, Symbol(colname)])
    #         @assert typ == cur_typ "Each modal frame with higher-than-zero dimension must have the same element types."
    #     end
    # end

    views = Vector{DataFrameRow}(undef, nrow(data))
    for i in 1:nrow(data)
        views[i] = view(data, i, :)
    end
    ModalFrame(dimension, data, views)
end

Base.size(mf::ModalFrame) = (nrow(mf.data),)
Base.getindex(mf::ModalFrame, i::Int) = mf.views[i]

struct ClassificationDataset
    ldim::Int
    hdim::Int
    frames::Vector{ModalFrame}
    instances::Vector{ModalInstance}
    classes::CategoricalArray
end

function ClassificationDataset(frames::Vector{ModalFrame}, classes::CategoricalArray)
    instances = ModalInstance[]

    n_instances = size(frames[1])

    for i in 1:length(frames)
        curr_n_instances = size(frames[i])
        @assert n_instances == curr_n_instances "Each frame must have the same number of (intraframe) instances!"
    end

    @assert n_instances[1] == length(classes) "Each instance must have one and only one class!"

    # for each instance
    for i in 1:n_instances[1]
        instance = ModalInstance()
        # for each frame
        for f in 1:length(frames)
            # push the instance
            push!(instance.rows, frames[f][i])
        end
        push!(instances, instance)
    end

    ldim, hdim = extrema([dimension(i) for i in frames])

    ClassificationDataset(ldim, hdim, frames, instances, classes)
end

instance(cs::ClassificationDataset, i::Int) = cs.instances[i] # TODO cs -> ds (either sup or unsupervised)
instance(cs::ClassificationDataset, i::Int, f::Int) = cs.instances[i][f]
function attributes(ds::ClassificationDataset)
    d = Dict{Int,Array{String,1}}()
    for (fid,frame) in enumerate(ds.frames)
        d[fid] = names(frame.data)
    end
    return d
end
attributes(ds::ClassificationDataset, fid::Int) = names(ds.frames[fid].data)

function Base.show(io::IO, ::MIME"text/plain", ds::ClassificationDataset)
    println(io, "Classification dataset with $(length(ds.instances)) instances")

    for (i, inst) in enumerate(ds.instances)
        println(io, "Instance #$(i):")
        println(io, instance(ds, i))
        println(io)
    end
end

# TODO
# possible test: all([auslan.instances[i].rows[1][attr] === auslan.frames[1].data[i,attr] for i in 1:length(auslan.instances) for attr in attributes(auslan, 1)])
function transform!(ds::ClassificationDataset, f::Function, fid::Int; kwargs...)
    for attr in attributes(ds, fid)
        for i in 1:length(nrow(ds.frames[fid].data))
            ds.frames[fid].data[i,attr] = f(ds.frames[fid].data[i,attr]; kwargs...)
        end
    end
end

function select_attributes!(ds::ClassificationDataset, frame_attributes::Dict{Int, Array{Int, 1}})
    frames = collect(keys(frame_attributes))

    for i in 1:length(frames)
        ds.frames[frames[i]].data = ds.frames[frames[i]].data[:, frame_attributes[frames[i]]]
    end
end

function resample(ds::ClassificationDataset, percentage::Real,
                seed::Int; balanced::Bool = true)
    n_instances = length(ds.instances)
    sample_size = n_instances * percentage / 100

    Random.seed!(seed)

    selected = Array([i for i in 1:n_instances])

    for i in 1:sample_size
        chosenLoc = rand(1:n_instances)
        chosen    = selected[chosenLoc]
        n_instances -= 1
        selected[chosenLoc] = selected[n_instances]
        selected[n_instances] = chosen
    end

    train = DataFrame()
    test  = DataFrame()
    train_classes = Int[]
    test_classes  = Int[]

    attributes = names(ds.frames[1].data)

    for attr in attributes
        insertcols!(train, attr => Array{Float64, 1}[])
        insertcols!(test, attr => Array{Float64, 1}[])
    end

    for i in 1:length(selected)
        if i > n_instances
            push!(train, ds.frames[1].data[selected[i], :])
            push!(train_classes, ds.classes[selected[i]])
        else
            push!(test, ds.frames[1].data[selected[i], :])
            push!(test_classes, ds.classes[selected[i]])
        end
    end

    train_ds = ClassificationDataset([ModalFrame(train)], CategoricalArray(train_classes))
    test_ds  = ClassificationDataset([ModalFrame(test)], CategoricalArray(test_classes))

    if balanced == true
        train_ds = down_sample(train_ds, seed)
        test_ds = down_sample(test_ds, seed)
    end
    return train_ds, test_ds
end

function down_sample(ds::ClassificationDataset, seed::Int)
    Random.seed!(seed)

    balanced_df = DataFrame()
    classes = Int[]

    attributes = names(ds.frames[1].data)
    for attr in attributes
        insertcols!(balanced_df, attr => Array{Float64, 1}[])
    end

    trip_ind = Int[]
    shut_ind = Int[]

    for i in 1:length(ds.classes)
        if ds.classes[i] == 1
            push!(trip_ind, i)
        else
            push!(shut_ind, i)
        end
    end

    trips = length(trip_ind)
    shutdowns = length(shut_ind)

    if trips < shutdowns
        selected = Array([i for i in 1:shutdowns])
        #rem = shutdowns

        for i in 1:trips
            chosenLoc = rand(1:shutdowns)
            chosen = shut_ind[chosenLoc]
            shut_ind[chosenLoc] = shut_ind[shutdowns]
            shutdowns-=1
            selected[i] = chosen
            #selected[chosenLoc] = selected[shutdowns]
            #selected[shutdowns] = chosen
        end

        for i in 1:trips
            push!(balanced_df, ds.frames[1].data[trip_ind[i], :])
            push!(classes, ds.classes[trip_ind[i]])
            push!(balanced_df, ds.frames[1].data[selected[i], :])
            push!(classes, ds.classes[selected[i]])
        end
    else
        selected = Array([i for i in 1:trips])

        for i in 1:shutdowns
            chosenLoc = rand(1:trips)
            chosen = shut_ind[chosenLoc]
            shut_ind[chosenLoc] = shut_ind[trips]
            trips-=1
            selected[i] = chosen
            #selected[chosenLoc] = selected[trips]
            #selected[trips] = chosen
        end

        for i in 1:shutdowns
            push!(balanced_df, ds.frames[1].data[shut_ind[i], :])
            push!(classes, ds.classes[shut_ind[i]])
            push!(balanced_df, ds.frames[1].data[selected[i], :])
            push!(classes, ds.classes[selected[i]])
        end
    end

    return ClassificationDataset([ModalFrame(balanced_df)], CategoricalArray(classes))
end

# function transform!(ds::ClassificationDataset, f::Function; kwargs...)

using MultivariateTimeSeries
X, y = read_data_labeled(joinpath(dirname(pathof(GBDTs)), "..", "data", "auslan_youtube8"));

# adesso voglia passare da MTS a ClassificationDataset

my_frame = DataFrame()

for attr in names(X)
    insertcols!(my_frame, attr => Array{Float64,1}[])
end

# for each instance
for i in 1:length(X)
    array_data = Array{Float64, 2}(X[i])
    push!(my_frame, [array_data[:, j] for j in 1:length(names(X))])
end

auslan = ClassificationDataset([ModalFrame(my_frame)], CategoricalArray(y))

# julia> p = MonteCarlo(30, 5)
# julia> using Random; Random.seed!(1)
# julia> model = induce_tree(grammar, :b, p, auslan.frames[1], y, 6);

"""

    if balanced == true
        #shutdowns = count(i -> i == 1, train_classes)
        #trips     = count(i -> i == 0, train_classes)

        ind_shut = Int[]
        ind_trip = Int[]
        #shutdowns = 1
        #trips = 1
        for i in 1:length(train_classes)
            if train_classes[i] == 1
                push!(ind_trip, i)
                #ind_shut[shutdowns] = i
                #shutdowns+=1
            else
                push!(ind_shut, i)
                #ind_trip[trips] = i
                #trips+=1
            end
        end

        shutdowns = length(ind_shut)
        trips = length(ind_trip)

        if shutdowns != trips
            if shutdowns < trips
                for i in 1:(trips-shutdowns)
                    chosen = rand(1:trips)
                    trips-=1
                    delete!(train, ind_trip[chosen])
                    splice!(train_classes, ind_trip[chosen])
                end
            else
                for i in 1:(shutdowns-trips)
                    chosen = rand(1:shutdowns)
                    shutdowns-=1
                    delete!(train, ind_shut[chosen])
                    splice!(train_classes, ind_shut[chosen])
                end
            end
        end

        #shutdowns = count(i -> i == 1, test_classes)
        #trips     = count(i -> i == 0, test_classes)
        #shutdowns = 1
        #trips = 1

        ind_shut = Int[]
        ind_trip = Int[]
        for i in 1:length(test_classes)
            if test_classes[i] == 1
                push!(ind_trip, i)
                #ind_shut[shutdowns] = i
                #shutdowns+=1
            else
                push!(ind_shut, i)
                #ind_trip[trips] = i
                #trips+=1
            end
        end

        if shutdowns != trips
            if shutdowns < trips
                for i in 1:(trips-shutdowns)
                    chosen = rand(1:trips)
                    trips-=1
                    delete!(test, ind_trip[chosen])
                    splice!(test_classes, ind_trip[chosen])
                end
            else
                for i in 1:(shutdowns-trips)
                    chosen = rand(1:shutdowns)
                    shutdowns-=1
                    delete!(train, ind_shut[chosen])
                    splice!(test_classes, ind_shut[chosen])
                end
            end
        end
    end
"""


text:
using Test
using FeatureEng

using Dates
using DataFrames

println("Running tests...")
println("Meep. Morp.")

@testset "Numeric" begin

    @testset "Scaling" begin
        data = [ 1:100; ]
        scale = fit_transform!(ScaleMinMax(),data)
        transformed = apply_transform(scale,data)
        @test minimum(transformed) == 0.0
        @test maximum(transformed) == 1.0

        data = [ 0, 5, 10 ]
        scale = fit_transform!(ScaleVariance(),data)
        transformed = apply_transform(scale,data)
        @test transformed == [-1.,0.,1.]

        data = [ 0, 0, 2^.5, 2^.5 ]
        scale = fit_transform!(ScaleL2(),data)
        transformed = apply_transform(scale,data)
        @test transformed == data ./ sum(data .^ 2) ^ .5

    end

    @testset "Transforming" begin
        @testset "Log Transform" begin
            data = ℯ .^ [ 1:10; ]
            @test isapprox(transformLog(data), [1:10;])

            for base = [ℯ,2,10]
                data = base .^ [ 1:10; ]
                @test isapprox(transformLog(data,base), [1:10;])
            end
        end

        @testset "Root Transform" begin
            for base = [ℯ,2,10]
                data = [1:10;] .^ base
                @test isapprox(transformRoot(data,base), [1:10;])
            end
        end

        @testset "BoxCox Transform" begin
            data = [ 1:10; ]

            @test isapprox(transformBoxCox(data,0.0), log.(data))
            @test isapprox(transformBoxCox(data,1.0), [ 0:9; ])

        end
    end

    @testset "Binning" begin
    end

    @testset "Interaction" begin
        df = DataFrame(
            x = 1:10,
            y = repeat([0,1],5)
        )
        col_names = ["x","x_x","x_y","y","y_y"]
        res = polynomial(df)
    end
end

@testset "Categorical" begin
    @testset "One-Hot Encoding" begin
        data = [3,1,2,4]
        @test encode_onehot(data) == DataFrame(Dict(
            "1" => [0,1,0,0],
            "2" => [0,0,1,0],
            "3" => [1,0,0,0],
            "4" => [0,0,0,1]
        ))
        @test encode_onehot(data,"col_") == DataFrame(Dict(
            "col_1" => [0,1,0,0],
            "col_2" => [0,0,1,0],
            "col_3" => [1,0,0,0],
            "col_4" => [0,0,0,1]
        ))
        @test encode_onehot(data,[1:6;]) == DataFrame(Dict(
            "1" => [0,1,0,0],
            "2" => [0,0,1,0],
            "3" => [1,0,0,0],
            "4" => [0,0,0,1],
            "5" => [0,0,0,0],
            "6" => [0,0,0,0]
        ))
        @test encode_onehot(data,[1:6;],"col_") == DataFrame(Dict(
            "col_1" => [0,1,0,0],
            "col_2" => [0,0,1,0],
            "col_3" => [1,0,0,0],
            "col_4" => [0,0,0,1],
            "col_5" => [0,0,0,0],
            "col_6" => [0,0,0,0],
        ))
        @test encode_onehot(["dog","cat","dog"]) == DataFrame(Dict(
            "cat" => [0,1,0],
            "dog" => [1,0,1]
        ))
    end

    @testset "Dummy Encoding" begin
    end

    @testset "Hash Encoding" begin
        data = [1:100:1000;]
        @test size(encode_hash(data,4)) == (10,4)
        @test encode_hash(data,4) == encode_hash(data,4)
    end
end

@testset "Datetime" begin
    @testset "Parse DateTime Strings" begin
        date_strings = [
            "2021-01-27 14:03:25",
            "1999-10-05 01:13:43",
            "abcdefg"
        ]
        datetimes = [
            DateTime(2021,1,27,14,03,25),
            DateTime(1999,10,5,1,13,43),
            missing
        ]
        @test all(a === b for (a,b) = zip(
            strp_datetimes(date_strings), datetimes))
    end
    @testset "Get Date Features" begin
        data = [
            DateTime(2021,1,27,14,3,25),
            DateTime(1999,10,5,1,13,43),
            DateTime(2010,6,12,11,0,0)
        ]
        @test extract_date_features(data) == DataFrame(
            year=[2021,1999,2010],
            month=["January","October","June"],
            dayofmonth=[27,5,12],
            dayofweek=["Wednesday","Tuesday","Saturday"],
            isweekend=[0,0,1],
            quarter=[1,4,2]
        )
    end
    @testset "Get Time Features" begin
        data = [
            DateTime(2021,1,27,14,3,25),
            DateTime(1999,10,5,1,13,43),
            DateTime(2010,6,12,11,0,0)
        ]
        @test extract_time_features(data) == DataFrame(
            hour=[14,1,11],
            minute=[3,13,0],
            second=[25.,43.,0.],
            isAM=[0,1,1]
        )
    end
    @testset "Get Date Features" begin
        data = [
            DateTime(2021,1,27,14,3,25),
            DateTime(1999,10,5,1,13,43),
            DateTime(2010,6,12,11,0,0)
        ]
        @test extract_datetime_features(data) == DataFrame(
            year=[2021,1999,2010],
            month=["January","October","June"],
            dayofmonth=[27,5,12],
            dayofweek=["Wednesday","Tuesday","Saturday"],
            isweekend=[0,0,1],
            quarter=[1,4,2],
            hour=[14,1,11],
            minute=[3,13,0],
            second=[25.,43.,0.],
            isAM=[0,1,1]
        )
    end
end


text:
include("src/H2O.jl")
using h2o

h2o.init()

h2o.clusterinfo()

h2o.get("Cloud")


text:
##
#
# PiBBP tests
#

using PiBBP
using Base.Test

#
# Particular digits
#

@test digpi(0)       == 2
@test digpi(1000000) == 6

#
# It actually computes fractional part of $\pi$.
#

@test pi - (3 + sum(digpi(n) / 16^(n+1) for n in 0:14)) ≈ 0
@test pi - (3 + sum(digpi(n) / big(16)^(n+1) for n in 0:64)) ≈ 0


text:
@testset "IO" begin

    msa_types = (
        Matrix{Residue},
        NamedResidueMatrix{Array{Residue,2}},
        MultipleSequenceAlignment,
        AnnotatedMultipleSequenceAlignment
        )

    # Sequence from PF09645
    F112_SSV1 = collect(".....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPD" *
                        "ECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ.....")

    # Test pfam stockholm parser using the 4 sequence full MSA for PF09645
    # > Order: Tree
    # > Inserts lower case
    # > Gaps as "." or "-" (mixed)
    # > Pfam version 28.0, based on UniProt release 2014_07
    @testset "Stockholm" begin

        pf09645_sto = joinpath(pwd(), "data", "PF09645_full.stockholm")

        @testset "Read" begin

            # TODO : @inferred read(pf09645_sto, Stockholm);

            @test isa(read(pf09645_sto, Stockholm), AnnotatedMultipleSequenceAlignment)
            for T in msa_types
                @test isa(read(pf09645_sto, Stockholm, T), T)
            end

            for T in msa_types
                @test read(pf09645_sto, Stockholm, T) == read(pf09645_sto, Stockholm)
            end
        end

        @testset "Output types" begin

            msa_objects = [ read(pf09645_sto, Stockholm, T)  for T in msa_types ]

            @testset "Sequence Names" begin

                default = ["1", "2", "3", "4"]
                pfnames = [ "C3N734_SULIY/1-95","H2C869_9CREN/7-104",
                            "Y070_ATV/2-70",    "F112_SSV1/3-112"   ]

                @test sequencenames(msa_objects[1]) == default
                for i in 2:4
                    @test sequencenames(msa_objects[i]) == pfnames
                end
            end

            @testset "Size" begin

                for msa in msa_objects
                    @test size(msa, 1) == 4
                    @test size(msa, 2) == sum(F112_SSV1 .!= Ref('.')) # without inserts
                    @test view(msa,4,:) == map(Residue, F112_SSV1[F112_SSV1 .!= Ref('.') ])
                end
            end
        end

        @testset "Annotations" begin

            msa = read(pf09645_sto, Stockholm)

            @test !isempty(msa)
            @test length(annotations(msa)) > 8 # 8 + modifications
            @test length(getannotcolumn(msa)) == 2
            @test length(getannotresidue(msa)) == 1
            @test getannotresidue(msa, "F112_SSV1/3-112", "SS") ==
                "X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
            @test getannotcolumn(msa, "seq_cons") ==
                "...NshphAclhaKILppKtElolEDIlAQFEISsosAYsI.+sL+hICEpH.-ECpsppKsRKTlhh.hKpEphppptpEp..ppItKIhsAp................"
            @test getannotsequence(msa, "F112_SSV1/3-112", "DR") == "PDB; 2VQC A; 4-73;"
        end

        @testset "String input" begin

            pfam_string = """
                # STOCKHOLM 1.0
                #=GS C3N734_SULIY/1-95   AC C3N734.1
                #=GS H2C869_9CREN/7-104  AC H2C869.1
                #=GS Y070_ATV/2-70       AC Q3V4T1.1
                #=GS F112_SSV1/3-112     AC P20220.1
                #=GS F112_SSV1/3-112     DR PDB; 2VQC A; 4-73;
                C3N734_SULIY/1-95              ...mp---NSYQMAEIMYKILQQKKEISLEDILAQFEISASTAYNVQRTLRMICEKHPDECEVQTKNRRTIFKWIKNEETTEEGQEE--QEIEKILNAQPAE-------------k....
                H2C869_9CREN/7-104             ...nk--LNDVQRAKLLVKILQAKGELDVYDIMLQFEISYTRAIPIMKLTRKICEAQ-EICTYDEKEHKLVSLNAKKEKVEQDEEENEREEIEKILDAH----------------trreq
                Y070_ATV/2-70                  qsvne-------VAQQLFSKLREKKEITAEDIIAIYNVTPSVAYAIFTVLKVMCQQHQGECQAIKRGRKTVI-------------------------------------------vskq.
                F112_SSV1/3-112                .....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ.....
                #=GR F112_SSV1/3-112     SS    .....X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX.....
                #=GC SS_cons                   .....X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX.....
                #=GC seq_cons                  ........NshphAclhaKILppKtElolEDIlAQFEISsosAYsI.+sL+hICEpH.-ECpsppKsRKTlhh.hKpEphppptpEp..ppItKIhsAp................h....
                //
                """

            @test parse(pfam_string, Stockholm) == read(pf09645_sto, Stockholm)

            msa = parse(pfam_string, Stockholm)
            @test !isempty(msa)
            @test length(annotations(msa)) > 8 # 8 + modifications
            @test length(getannotcolumn(msa)) == 2
            @test length(getannotresidue(msa)) == 1
        end

        @testset "File write Matrix{Residue}" begin

            path = tempdir()
            tmp_file = joinpath(path, ".tmp.stockholm")
            try
                msa = read(pf09645_sto, Stockholm, Matrix{Residue})
                write(tmp_file, msa, Stockholm)
                @test read(tmp_file, Stockholm, Matrix{Residue}) == msa
            finally
                if isfile(tmp_file)
                    rm(tmp_file)
                end
            end
        end

        @testset "Keep insert columns" begin

            msa = read(pf09645_sto, Stockholm, keepinserts=true)
            # Aligned columns
            @test (collect(getannotcolumn(msa, "Aligned")) .== Ref('1')) == (F112_SSV1 .!= Ref('.'))
            @test stringsequence(msa, 1) == replace(uppercase(
                "...mp---NSYQMAEIMYKILQQKKEISLEDILAQFEISASTAYNVQRTLRMICEKHPDECEVQTKNRRTIFKWIKNEETTEEGQEE--QEIEKILNAQPAE-------------k...."
                ), '.' => '-')
        end
    end

    @testset "FASTA" begin

        pf09645_fas = joinpath(pwd(), "data", "PF09645_full.fasta.gz")
        gaoetal2011 = joinpath(pwd(), "data", "Gaoetal2011.fasta")

        @testset "Read" begin

            @test isa(read(pf09645_fas, FASTA), AnnotatedMultipleSequenceAlignment)
            for T in msa_types
                @test isa(read(pf09645_fas, FASTA, T), T)
            end

            for T in msa_types
                @test read(pf09645_fas, FASTA, T) == read(pf09645_fas, FASTA)
            end

            @testset "Download" begin

                @test read(gaoetal2011, FASTA) ==
                    read("https://raw.githubusercontent.com/diegozea/MIToS.jl/master/test/data/Gaoetal2011.fasta",
                        FASTA)
                @test read(pf09645_fas, FASTA) ==
                    read("https://raw.githubusercontent.com/diegozea/MIToS.jl/master/test/data/PF09645_full.fasta.gz",
                        FASTA)
            end
        end

        @testset "Output types" begin

            gaoetal_msas = [ read(gaoetal2011, FASTA, T)  for T in msa_types ]
            pfam_msas  = [ read(pf09645_fas, FASTA, T)  for T in msa_types ]

            @testset "Sequence Names" begin

                @test sequencenames(gaoetal_msas[1]) == ["1", "2", "3", "4", "5", "6"]
                @test sequencenames(pfam_msas[1])  == ["1", "2", "3", "4"]
                for i in 2:4
                    @test sequencenames(gaoetal_msas[i]) ==
                        ["SEQ1", "SEQ2", "SEQ3", "SEQ4", "SEQ5", "SEQ6"]
                    @test sequencenames(pfam_msas[i]) == [  "C3N734_SULIY/1-95",
                                                            "H2C869_9CREN/7-104",
                                                            "Y070_ATV/2-70",
                                                            "F112_SSV1/3-112"   ]
                end
            end

            @testset "Residues" begin

                list_seq = [res"DAWAEE",
                            res"DAWAEF",
                            res"DAWAED",
                            res"DAYCMD",
                            res"DAYCMT",
                            res"DAYCMT" ]
                residues = permutedims(hcat(list_seq...), [2,1])

                for msa in gaoetal_msas

                    @test msa == residues
                    @test isa(getresidues(msa), Matrix{Residue})
                    @test getresidues(msa) == residues
                    @test getresiduesequences(msa) == list_seq
                    @test stringsequence(msa, 1) == "DAWAEE"
                end

                for msa in gaoetal_msas[2:end]
                    # getsequence returns a matrix, use vec(seq) or dropdims(seq,  dims=1)
                    # to get a vector:
                    @test vec(getresidues(getsequence(msa, 1))) == res"DAWAEE"
                    @test dropdims(getresidues(getsequence(msa, 1)),  dims=1) == res"DAWAEE"
                end
            end

            @testset "Size" begin

                for msa in gaoetal_msas
                    @test size(msa, 1) == 6
                    @test size(msa, 2) == 6
                    @test view(msa,4,:) == res"DAYCMD"
                end

                for msa in pfam_msas
                    @test size(msa, 1) == 4
                    @test size(msa, 2) == sum(F112_SSV1 .!= Ref('.')) # without inserts
                    @test view(msa,4,:) == map(Residue, F112_SSV1[ F112_SSV1 .!= Ref('.') ])
                end
            end
        end

        @testset "String input/output" begin

            msa = read(gaoetal2011, FASTA)
            fasta_string = """
                >SEQ1
                DAWAEE
                >SEQ2
                DAWAEF
                >SEQ3
                DAWAED
                >SEQ4
                DAYCMD
                >SEQ5
                DAYCMT
                >SEQ6
                DAYCMT
                """

            @test parse(fasta_string, FASTA) == msa

            out = IOBuffer()
            print(out, msa, FASTA)
            @test String(take!(out)) == fasta_string
        end

        @testset "File input/output" begin

            msa = read(gaoetal2011, FASTA)

            path = tempdir()
            uncompressed = joinpath(path, ".tmp.fasta")
            try
                write(uncompressed, msa, FASTA)
                @test read(uncompressed, FASTA) == msa
            finally
                if isfile(uncompressed)
                    rm(uncompressed)
                end
            end
            compressed = joinpath(path, ".tmp.fasta.gz")
            try
                write(compressed, msa, FASTA)
                @test read(compressed, FASTA) == msa
            finally
                if isfile(compressed)
                    rm(compressed)
                end
            end
        end

        @testset "Non standard residues and mapping" begin

            seqs = read(joinpath(pwd(), "data", "alphabet.fasta"), FASTA,
                        generatemapping=true)

            @test vec(seqs[1,:]) == res"ARNDCQEGHILKMFPSTWYV"
            for i in 2:nsequences(seqs)
                @test vec(seqs[i,:]) == res"ARXDCQEGHILKMFPSTWYV"
                @test getsequencemapping(seqs,1) == getsequencemapping(seqs,i)
            end
        end
    end

    @testset "Raw" begin

        # AnnotatedMultipleSequenceAlignment
        raw = read(joinpath(pwd(), "data", "gaps.txt"), Raw)
        mat = getresidues(raw)
        raw_string = """THAYQAIHQV
                        THAYQAIHQ-
                        THAYQAIH--
                        THAYQAI---
                        THAYQA----
                        THAYQ-----
                        THAY------
                        THA-------
                        TH--------
                        T---------
                        """
        @testset "Parse" begin

            @test stringsequence(raw, "1") == "THAYQAIHQV"
            @test stringsequence(raw, 10)  == "T---------"

            for i in 1:10
                @test stringsequence(mat, i)  == stringsequence(raw, i)
            end

            @test parse(raw_string, Raw) == raw
        end

        @testset "Print" begin

            buffer = IOBuffer()
            print(buffer, raw, Raw)
            @test String(take!(buffer)) == raw_string

            print(buffer, mat, Raw)
            @test String(take!(buffer)) == raw_string
        end

        @testset "Stats" begin

            @test gapfraction(mat) == 0.45
            @test vec(gapfraction(mat, 1)) ≈ [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]
            @test vec(gapfraction(mat, 2)) ≈ [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]

            @test residuefraction(mat) == 0.55
            @test vec(residuefraction(mat, 1)) ≈ [1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1]
            @test vec(residuefraction(mat, 2)) ≈ [1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1]

            @test vec(coverage(mat)) ≈ [1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1]

            @test vec(columngapfraction(mat)) ≈ [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]
        end

        @testset "Reference and gapstrip" begin

            gs = gapstrip(mat, coveragelimit=0.5, gaplimit=0.5)
            @test vec(getsequence(gs, 1)) == res"THAYQAIH"
            @test ncolumns(gs) == 8
            @test nsequences(gs) == 6

            ref = setreference!(copy(mat), 2)
            @test vec(getsequence(ref, 1)) == res"THAYQAIHQ-"
            ref = adjustreference(ref)
            @test vec(getsequence(ref, 1)) == res"THAYQAIHQ"
        end
    end

    @testset "NBRF/PIR" begin

        # Example NBRF file: http://iubio.bio.indiana.edu/soft/molbio/readseq/classic/src/Formats
        # "The sequence is free format and may be interrupted by blanks for ease of reading"
        example = joinpath(pwd(), "data", "example.nbrf")
        # Alignment file (PIR) from https://salilab.org/modeller/9v7/manual/node445.html
        modeller = joinpath(pwd(), "data", "modeller.pir.gz")
        # http://emboss.sourceforge.net/docs/themes/seqformats/NbrfFormat.html
        # "sequence may contain punctuation symbols to indicate various degrees of
        # reliability of the data"
        emboss = joinpath(pwd(), "data", "emboss.pir")
        # http://caps.ncbs.res.in/pass2v3/pir.html
        # Example from pass2 with spaces added at the end of the id lines.
        pass2 = joinpath(pwd(), "data", "pass2.pir")

        @testset "Read" begin

            @test isa(read(modeller, PIR), AnnotatedMultipleSequenceAlignment)
            for T in msa_types
                @test isa(read(modeller, PIR, T), T)
            end

            msa = read(modeller, PIR)
            @test size(msa) == (2, 106)
        end

        @testset "Print" begin

            msa = read(example, PIR)

            @test stringsequence(msa, 1) == "MTNIRKSHPLFKIINHSFIDLPAPSVTHICRDVNYGWLIRYTWIGGQPVEHPFIIIGQLASISYFSIILILMPISGIVEDKMLKWN"

            out = IOBuffer()
            print(out, msa, PIR)
            @test String(take!(out)) == """
                >P1;CBRT
                Cytochrome b - Rat mitochondrion (SGC1)
                MTNIRKSHPLFKIINHSFIDLPAPSVTHICRDVNYGWLIRYTWIGGQPVEHPFIIIGQLASISYFSIILILMPISGIVED
                KMLKWN*
                """

            out = IOBuffer()
            print(out, msa.matrix, PIR) # NamedResidueMatrix{Array{Residue,2}}
            @test String(take!(out)) == """
                >XX;CBRT

                MTNIRKSHPLFKIINHSFIDLPAPSVTHICRDVNYGWLIRYTWIGGQPVEHPFIIIGQLASISYFSIILILMPISGIVED
                KMLKWN*
                """

            out = IOBuffer()
            print(out, msa.matrix.array, PIR) # Matrix{Residue}
            @test String(take!(out)) == """
                >XX;1

                MTNIRKSHPLFKIINHSFIDLPAPSVTHICRDVNYGWLIRYTWIGGQPVEHPFIIIGQLASISYFSIILILMPISGIVED
                KMLKWN*
                """

        end

        @testset "Gaps" begin

            msa = read(modeller, PIR)
            @test gapfraction(msa, 2) ≈ [0.0, 0.49056603773584906]
            @test getannotsequence(msa, "5fd1", "Type") == "P1"
            @test getannotsequence(msa, "5fd1",
                "Title") == "structureX:5fd1:1    :A:106  :A:ferredoxin:Azotobacter vinelandii: 1.90: 0.19"
        end

        @testset "String input/output" begin

            emboss_string = """
                >P1;AZBR
                finger protein zfpA - turnip fern chloroplast
                GDVE(G.K.G.I.F=T,M,C.S.Q,C.H.V,E.K.G.G.K.H)
                FTGPNLHGLFGRK.TGQAVGYSYTAANK.NK.GIIWGDDTLM
                EYLENPK.RYIPGTK.MVFTGLSK.YRE
                RTNLIAYLK.EK.TAA*
                """
            # MIToS always reads . as -
            msa = read(emboss, PIR, deletefullgaps=false)
            @test parse(emboss_string, PIR, deletefullgaps=false) == msa

            out = IOBuffer()
            print(out, msa, PIR)
            @test String(take!(out)) == """
                >P1;AZBR
                finger protein zfpA - turnip fern chloroplast
                GDVEG-K-G-I-FTMC-S-QC-H-VE-K-G-G-K-HFTGPNLHGLFGRK-TGQAVGYSYTAANK-NK-GIIWGDDTLMEY
                LENPK-RYIPGTK-MVFTGLSK-YRERTNLIAYLK-EK-TAA*
                """
        end

        @testset "Spaces at the end of the id line" begin
            lines = readlines(pass2)
            @test lines[1] == ">P1;1bbha- "
            @test lines[5] == ">P1;1cpq-- "
            @test lines[9] == ">P1;256bb- "
            msa = read(pass2, PIR)
            @test sequencenames(msa) == ["1bbha-", "1cpq--", "256bb-"]
        end

        @testset "Duplicated identifiers" begin

            duplicated_ids_file = joinpath(pwd(), "data", "duplicated_ids.pir")
            @test_throws ArgumentError read(duplicated_ids_file, PIR)
        end
    end
end


text:
"""

Compare two vectors of elements using Base.isapprox for each element's
ray matrix (ABCD entries). Does consequently not consider one
discretization of element FreeSpace different from another, or one
realization of an imaging system from another as long as both achieve
(within tolerances) the same imaging.

!!! note
    The `atol` (absolute tolerance) parameter can be used but is
    typically nonsensical as it will be used for each of the
    ray matrix entries ABCD which usually differ vastly in magnitude.

"""
Base.isapprox(
    a::Vector{<:Element}, b::Vector{<:Element}; kwargs...
) = isapprox(RTM(a), RTM(b); kwargs...)


text:
function findFirstSmaller(x::Float64,v::Array{Float64,1})

    l = length(v)
    if l == 129         # 8bit case
        idx = 64        # starting index (l-1)/2
        n = 6           #
    elseif l == 32769   # 16bit case
        idx = 16384     # starting index (l-1)/2
        n = 14
    else
        throw(BoundsError())
    end

    if x > v[end-1]     # floatmax and infinity case
        return l
    end

    # binary tree search
    for i in 1:n
        if x < v[idx]
            idx -= 2^(n-i)
        else
            idx += 2^(n-i)
        end
    end

    # split off the i = n+1 case
    if x >= v[idx]
        idx += 1
    end

    return idx
end