text:
"""
REINFORCE baseline
------------------

REINFORCE with the state values V as a baseline

See REINFORCE.jl for some more implementational info.
See learner/REINFORCEbaselinelearner.jl for the core implementation.

Note
----
get_all_puzzles(), get_all_episodes_in_puzzle(), get_data_by_puzzle_episode()
are in rl_data_toolsforalgos.jl

"""

function getLL_REINFORCE_baseline(data::Array{Int,2}, alpha::Float64,
                        gamma::Float64, temp::Float64, alpha_baseline::Float64;
                        do_debug_print = false)
    data_stats, nr_states, nr_actions, nr_samples = get_data_stats_and_dimensions(data)
    if data_stats.nr_subjects != 1
        error("Invalid argument: getLL_REINFORCE_baseline supports only one subject")
    end
    if nr_actions!=2
        error("getLL_REINFORCE_baseline supports two actions only")
    end
    all_puzzles = get_all_puzzles(data)

    LL_unbiased = 0.
    probs_unbiased = zeros(Float64, nr_samples)
    TD_error_record = zeros(Float64, nr_samples)

    Δpolicy_record = zeros(Float64, nr_samples)
    PolicyParametersImbalance_record = zeros(Float64, nr_samples)
    PolicyImbalance_record = zeros(Float64, nr_samples)
    PolicyParametersImbalanceAllstates_record = Array{Array{Float64,1}, 1}(undef, nr_samples)
    [PolicyParametersImbalanceAllstates_record[i] = zeros(Float64, nr_states) for i in 1:nr_samples]
    PolicyParametersDiffChosenVsUnchosen_record = zeros(Float64, nr_samples)
    PolicyDiffChosenVsUnchosen_record = zeros(Float64, nr_samples)

    modulatingfactor = 0.
    learner = []
    action_counter = 0
    for i_puzzle in all_puzzles
        # Puzzle = independent block of the experiment.
        # The code was made for general purposes, in this experiment we
        # have only ONE puzzle.
        learner = REINFORCEbaseline(nr_states = nr_states,
                        nr_actions = nr_actions,
                        α = alpha,
                        γ = gamma,
                        α_baseline = alpha_baseline,
                        do_debug_print = do_debug_print)

        all_epi_in_puzzle = get_all_episodes_in_puzzle(data, i_puzzle)
        for i_episode in all_epi_in_puzzle
            episode_data = get_data_by_puzzle_episode(data, i_puzzle, i_episode)
            path_length = size(episode_data)[1]  # = T
            if do_debug_print
                println("i_puzzle:$i_puzzle | i_episode:$i_episode | path_length:$path_length")
            end
            reward = 0. # we have a single reward at goal state. If the episode was
            # not completed, this code still correctly works (updates are 0)
            for t in 1:path_length
                action_counter += 1
                state = episode_data[t, Int(c_state)]
                action = episode_data[t, Int(c_action)]
                reward = episode_data[t, Int(c_reward)]
                next_state = episode_data[t, Int(c_nextState)]
                done = reward>0 ? true : false

                p_actions = getSoftmaxProbs([learner.θ[state,1],learner.θ[state,2]], temp, [0., 0.])

                LL_unbiased += log(p_actions[action])
                probs_unbiased[action_counter] = copy(p_actions[action])

                PolicyParametersImbalance_record[action_counter] = abs(learner.θ[state, 1] - learner.θ[state, 2])
                PolicyImbalance_record[action_counter] = abs(p_actions[1] - p_actions[2])
                PolicyParametersImbalanceAllstates_record[action_counter] = abs.(learner.θ[:, 1] - learner.θ[:, 2])
                PolicyParametersDiffChosenVsUnchosen_record[action_counter] = learner.θ[state, action] - learner.θ[state, mod(action,2)+1]
                PolicyDiffChosenVsUnchosen_record[action_counter] = p_actions[action] - p_actions[mod(action,2)+1]

                θ_before = deepcopy(learner.θ)

                update!(learner, state, action, reward, p_actions[action], modulatingfactor)
                TD_error_record[action_counter] = copy(learner.target[1])
                Δpolicy_record[action_counter] = sum(abs.(learner.θ .- θ_before))

                if do_debug_print
                    println("t:$t | s:$state | a:$action | R:$reward | s':$next_state | p(s,a):$(round(p_action,digits=2))")
                end
            end # path_length (in single episode)
        end # episode
    end # puzzzle

  LL = LL_unbiased
  trial_by_trial_p_sa = (probs_unbiased)

  algo = "REINFORCE with baseline"
  param_v = [alpha, gamma, temp, alpha_baseline]
  param_str = "alpha=$(round(alpha,digits=4)), gamma=$(round(gamma,digits=4)),
                temp=$(round(temp,digits=4)), alpha_baseline=$(round(alpha_baseline,digits=4))"

  results_dict = Dict{String,Any}("final_Theta" => learner.θ,
                                    "DeltaPolicy" => Δpolicy_record,
                                    "PolicyParametersImbalance" => PolicyParametersImbalance_record,
                                    "PolicyImbalance" => PolicyImbalance_record,
                                    "PolicyParametersImbalanceAllstates" => PolicyParametersImbalanceAllstates_record,
                                    "PolicyParametersDiffChosenVsUnchosen" => PolicyParametersDiffChosenVsUnchosen_record,
                                    "PolicyDiffChosenVsUnchosen" => PolicyDiffChosenVsUnchosen_record)

  signals = Step_by_step_signals(
    data_stats.all_subject_ids[1],
    TD_error_record,
    [], # no SPE
    trial_by_trial_p_sa,
    LL,
    data_stats,
    algo,
    param_v,
    param_str,
    99, #biased_action_idx, dummy value
    results_dict)
  return LL, signals

end


text:
#=
One-dimensional Kriging method, following this paper:
"A Taxonomy of Global Optimization Methods Based on Response Surfaces"
by DONALD R. JONES
=#
mutable struct Kriging{X,Y,L,U,P,T,M,B,S,R} <: AbstractSurrogate
    x::X
    y::Y
    lb::L
    ub::U
    p::P
    theta::T
    mu::M
    b::B
    sigma::S
    inverse_of_R::R
 end

 """
 Gives the current estimate for array 'val' with respect to the Kriging object k.
 """
 function (k::Kriging)(val)
     d = length(val)
     n = length(k.x)
     return k.mu + sum(k.b[i] * exp( -sum(k.theta[j]*norm(val[j]-k.x[i][j])^k.p[j] for j = 1:d) ) for i = 1:n)
 end

 """
     Returns sqrt of expected mean_squared_error error at the point.
 """
 function std_error_at_point(k::Kriging,val)
     n = length(k.x)
     d = length(k.x[1])
     r = zeros(eltype(k.x[1]),n,1)
     @inbounds for i = 1:n
         sum = zero(eltype(k.x[1]))
         for l = 1:d
             sum = sum + k.theta[l]*norm(val[l]-k.x[i][l])^(k.p[l])
         end
         r[i] = exp(-sum)
     end
     one = ones(eltype(k.x[1]),n,1)
     one_t = one'
     a = r'*k.inverse_of_R*r
     a = a[1]
     b = one_t*k.inverse_of_R*one
     b = b[1]
     mean_squared_error = k.sigma*(1 - a + (1-a)^2/(b))
     return sqrt(abs(mean_squared_error))
 end

 """
 Gives the current estimate for 'val' with respect to the Kriging object k.
 """
 function (k::Kriging)(val::Number)
     phi = z -> exp(-(abs(z))^k.p)
     n = length(k.x)
     prediction = zero(eltype(k.x[1]))
     for i = 1:n
         prediction = prediction + k.b[i]*phi(val-k.x[i])
     end
     prediction = k.mu + prediction
     return prediction
 end

"""
    Returns sqrt of expected mean_squared_error error at the point.
"""
function std_error_at_point(k::Kriging,val::Number)
    phi(z) = exp(-(abs(z))^k.p)
    n = length(k.x)
    r = zeros(eltype(k.x[1]),n,1)
    @inbounds for i = 1:n
        r[i] = phi(val - k.x[i])
    end
    one = ones(eltype(k.x[1]),n,1)
    one_t = one'
    a = r'*k.inverse_of_R*r
    a = a[1]
    b = one_t*k.inverse_of_R*one
    b = b[1]
    mean_squared_error = k.sigma*(1 - a + (1-a)^2/(b))
    return sqrt(abs(mean_squared_error))
end

"""
    Kriging(x,y,lb::Number,ub::Number;p::Number=1.0)

Constructor for type Kriging.

#Arguments:
-(x,y): sampled points
-p: value between 0 and 2 modelling the
   smoothness of the function being approximated, 0-> rough  2-> C^infinity
"""
function Kriging(x,y,lb::Number,ub::Number;p=1.0)
    if length(x) != length(unique(x))
        println("There exists a repetion in the samples, cannot build Kriging.")
        return
    end
    mu,b,sigma,inverse_of_R = _calc_kriging_coeffs(x,y,p)
    theta = 1.0
    Kriging(x,y,lb,ub,p,theta,mu,b,sigma,inverse_of_R)
end

function _calc_kriging_coeffs(x,y,p::Number)
    n = length(x)
    R = zeros(eltype(x[1]), n, n)
    @inbounds for i = 1:n
        for j = 1:n
            R[i,j] = exp(-abs(x[i]-x[j])^p)
        end
    end
    one = ones(eltype(x[1]),n,1)
    one_t = one'
    inverse_of_R = inv(R)
    mu = (one_t*inverse_of_R*y)/(one_t*inverse_of_R*one)
    b = inverse_of_R*(y-one*mu)
    sigma = ((y-one*mu)' * inverse_of_R * (y - one*mu))/n
    mu[1], b, sigma[1],inverse_of_R
end

"""
    Kriging(x,y,lb,ub;p=collect(one.(x[1])),theta=collect(one.(x[1])))

Constructor for Kriging surrogate.

- (x,y): sampled points
- p: array of values 0<=p<2 modeling the
     smoothness of the function being approximated in the i-th variable.
     low p -> rough, high p -> smooth
- theta: array of values > 0 modeling how much the function is
          changing in the i-th variable.
"""
function Kriging(x,y,lb,ub;p=collect(one.(x[1])),theta=collect(one.(x[1])))
    if length(x) != length(unique(x))
        println("There exists a repetition in the samples, cannot build Kriging.")
        return
    end
    mu,b,sigma,inverse_of_R = _calc_kriging_coeffs(x,y,p,theta)
    Kriging(x,y,lb,ub,p,theta,mu,b,sigma,inverse_of_R)
end

function _calc_kriging_coeffs(x,y,p,theta)
    n = length(x)
    d = length(x[1])
    R = zeros(float(eltype(x[1])), n, n)
    @inbounds for i = 1:n
        for j = 1:n
            sum = zero(eltype(x[1]))
            for l = 1:d
            sum = sum + theta[l]*norm(x[i][l]-x[j][l])^p[l]
            end
            R[i,j] = exp(-sum)
        end
    end
    one = ones(n,1)
    one_t = one'
    inverse_of_R = inv(R)
    mu = (one_t*inverse_of_R*y)/(one_t*inverse_of_R*one)
    b = inverse_of_R*(y-one*mu)
    sigma = ((y-one*mu)' * inverse_of_R * (y - one*mu))/n
    mu[1], b, sigma[1],inverse_of_R
end

"""
    add_point!(k::Kriging,new_x,new_y)

Adds the new point and its respective value to the sample points.
Warning: If you are just adding a single point, you have to wrap it with [].
Returns the updated Kriging model.

"""
function add_point!(k::Kriging,new_x,new_y)
    if new_x in k.x
        println("Adding a sample that already exists, cannot build Kriging.")
        return
    end
    if (length(new_x) == 1 && length(new_x[1]) == 1) || ( length(new_x) > 1 && length(new_x[1]) == 1 && length(k.theta)>1)
        push!(k.x,new_x)
        push!(k.y,new_y)
    else
        append!(k.x,new_x)
        append!(k.y,new_y)
    end
    k.mu,k.b,k.sigma,k.inverse_of_R = _calc_kriging_coeffs(k.x,k.y,k.p,k.theta)
    nothing
end

ext: jl


text:
module PhysicalBasesTests

using FillArrays
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Arrays
using Gridap.Fields
using Gridap.FESpaces
using Gridap.Polynomials
using Gridap.Integration
using Gridap.CellData
using Gridap.TensorValues
using Test

# Start with a PhysicalSpaceCellBasis
a = 1
b = ( a == 1 )

# domain = (0,1)
# partition = (3,)
domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
order = 1
# order = 2

trian = get_triangulation(model)
quad = CellQuadrature(trian,order)
q = get_coordinates(quad)

polytopes = get_polytopes(model)
cell_to_ctype = get_cell_type(model)

grid = get_grid(model)
cell_map = get_cell_map(grid)

# Test against the ref approach...

T = Float64
reffes = [LagrangianRefFE(T,p,order) for p in polytopes]

dof_bases = map(get_dof_basis,reffes)

cell_dof_basis = FESpaces._cell_dof_basis_physical_space(dof_bases,cell_to_ctype,cell_map)
cell_dof_basis = GenericCellDofBasis(Val{false}(),cell_dof_basis)

prebasis =  map(get_prebasis,reffes)
cell_prebasis = CompressedArray(prebasis,cell_to_ctype)

ncells = num_cells(model)
ndofs = num_dofs(first(reffes))
cell_prebasis_new = GenericCellField(cell_prebasis,cell_map,Val{false}(),Fill((Base.OneTo(ndofs),),ncells),Val((:,)))

# cell_matrix = evaluate(cell_dof_basis,cell_prebasis)

# cell_shapefuns = _cell_shape_functions_physical_space(cell_prebasis,cell_dof_basis,cell_map)

psfs, x  = FESpaces.compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
sfs, x  = FESpaces.compute_cell_space(reffes, cell_to_ctype, cell_map)

# T = VectorValue{2,Float64}
# reffes = [LagrangianRefFE(T,p,order) for p in polytopes]

r = evaluate(sfs,q)
rg = evaluate(gradient(sfs),q)
rp = evaluate(psfs,q)
rgp = evaluate(gradient(psfs),q)

@test all([ rg[i] ≈ rgp[i] for i in 1:length(rg) ])
@test all([ r[i] ≈ rp[i] for i in 1:length(rg) ])

T = VectorValue{2,Float64}
reffes = [LagrangianRefFE(T,p,order) for p in polytopes]

psfs, x  = FESpaces.compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
sfs, x  = FESpaces.compute_cell_space(reffes, cell_to_ctype, cell_map)

r = evaluate(sfs,q)
rg = evaluate(gradient(sfs),q)
rp = evaluate(psfs,q)
rgp = evaluate(gradient(psfs),q)

@test all([ rg[i] ≈ rgp[i] for i in 1:length(rg) ])
@test all([ r[i] ≈ rp[i] for i in 1:length(rg) ])

func(x) = x

cell_field = convert_to_cell_field(func,cell_map)
isa(cell_field,CellField)

# Now RT elements

T = Float64
order = 0
reffes = [RaviartThomasRefFE(T,p,order) for p in polytopes]

psfs, dofp  = FESpaces.compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
sfs, dof  = FESpaces.compute_cell_space(reffes, cell_to_ctype, cell_map)

r = evaluate(sfs,q)
rg = evaluate(gradient(sfs),q)
rp = evaluate(psfs,q)
rgp = evaluate(gradient(psfs),q)

@test all([ r[i] ≈ rp[i] for i in 1:length(rp) ])
@test all([ rg[i] ≈ rgp[i] for i in 1:length(rg) ])

end #module

# # If I want new evaluation...
# function kernel_evaluate(k::typeof{change_basis},x,cell_prebasis,cell_matrix_inv)
#    cell_prebasis_x = evaluate_field_array(cell_prebasis,x)
#   apply(mul,cell_prebasis_x,cell_prebasis,cell_matrix_inv)
# end
# function apply_gradient(k::typeof(change_basis),cell_prebasis,cell_matrix_inv)
#    cell_prebasis_grad = gradient(cell_prebasis)
#    apply(change_basis,cell_prebasis_grad,cell_matrix_inv)
# end
# ##
# # Optimisation : evaluate_field_array for AbstractArray with FieldLike
# # Define a new kernel that better treats the inverse
# struct InvKernel <: Kernel end
# function kernel_cache(k::InvKernel,mat)
# end
# function apply_kernel!(cache,k::InvKernel,mat)
# end
# function kernel_cache(k::InvKernel,mat)
# CachedArray(copy(mat))
# end
# function apply_kernel!(cache,k::InvKernel,mat)
#   setsize!(cache,size(mat))
#   m = cache.array
#   fill!(m,zero(m))
#   for i:size(m,1); m[i] = 1; end
#   ldiv!(mat,m)
#   m
# end
# k = InvKernel()
#
# isa(cell_prebasis,CellBasis)
#
# change_basis(cell_prebasis[1],cell_matrix_inv[1])


text:
using ITensors
using ITensorNetworkAD
using OptimKit
using Random
using Zygote

Random.seed!(1243)

N = 6
s = siteinds("S=1/2", N; conserve_qns=true)
os = OpSum()
for n in 1:(N - 1)
  os .+= 0.5, "S+", n, "S-", n + 1
  os .+= 0.5, "S-", n, "S+", n + 1
  os .+= "Sz", n, "Sz", n + 1
end
Hmpo = MPO(os, s)
ψ₀mps = randomMPS(s, n -> isodd(n) ? "↑" : "↓")
H = prod(Hmpo)
ψ₀ = prod(ψ₀mps)

# The Rayleigh quotient to minimize
function E(H::ITensor, ψ::ITensor)
  ψdag = dag(ψ)
  return (ψdag' * H * ψ)[] / (ψdag * ψ)[]
end

@show flux(ψ₀)
@show E(H, ψ₀)

E(ψ::ITensor) = E(H, ψ)
∇E(ψ::ITensor) = E'(ψ)
fg(ψ::ITensor) = (E(ψ), ∇E(ψ))

linesearch = HagerZhangLineSearch(;
  c₁=0.1, c₂=0.9, ϵ=1e-6, θ=1 / 2, γ=2 / 3, ρ=5.0, verbosity=0
)
algorithm = LBFGS(3; maxiter=30, gradtol=1e-8, linesearch=linesearch)
ψ, fψ, gψ, numfg, normgradhistory = optimize(fg, ψ₀, algorithm)

display(normgradhistory)

@show E(H, ψ)

D, _ = eigen(H; ishermitian=true)
@show minimum(D)


text:
export RK3

"""
    RK3()

The classical [Runge-Kutta](https://en.wikipedia.org/wiki/Runge–Kutta_methods)
predictor of order 3.
"""
struct RK3 <: AbstractPredictor end
struct RK3Cache{T} <: AbstractPredictorCache
    dt::Vector{T}
    k1::Vector{T}
    k2::Vector{T}
    k3::Vector{T}
    k4::Vector{T}
end

function cache(::RK3, H, x, ẋ, t)
    RK3Cache(dt(H, x, t), copy(ẋ), copy(ẋ), copy(ẋ), copy(ẋ))
end
#
function predict!(
    xnext,
    cache::RK3Cache,
    H::HomotopyWithCache,
    x,
    t,
    Δt,
    ẋ,
    Jac::JacobianMonitor,
)
    dt, mk₁, mk₂, mk₃ = cache.dt, cache.k1, cache.k2, cache.k3
    n = length(xnext)
    @inbounds for i = 1:n
        mk₁[i] = -ẋ[i]
        xnext[i] = x[i] - 0.5Δt * mk₁[i]
    end
    minus_ẋ!(mk₂, H, xnext, t + 0.5Δt, Jac, dt)

    @inbounds for i = 1:n
        xnext[i] = x[i] + Δt * mk₁[i] - 2Δt * mk₂[i]
    end
    minus_ẋ!(mk₃, H, xnext, t + Δt, Jac, dt)
    @inbounds for i = 1:n
        xnext[i] = x[i] - Δt * mk₃[i]
    end

    @inbounds for i = 1:n
        xnext[i] = x[i] - Δt / 6 * (mk₁[i] + 4 * mk₂[i] + mk₃[i])
    end

    nothing
end

order(::RK3Cache) = 4


text:
@testset "Access and Iterations" begin
    dna_kmer = DNAKmer("ACTG")
    rna_kmer = RNAKmer("ACUG")

    @testset "Access DNA Kmer" begin
        @test dna_kmer[1] == DNA_A
        @test dna_kmer[2] == DNA_C
        @test dna_kmer[3] == DNA_T
        @test dna_kmer[4] == DNA_G

        # Access indexes out of bounds
        @test_throws Exception dna_kmer[-1]
        @test_throws Exception dna_kmer[0]
        @test_throws Exception dna_kmer[5]
        @test_throws Exception getindex(dna_kmer,-1)
        @test_throws Exception getindex(dna_kmer, 0)
        @test_throws Exception getindex(dna_kmer, 5)
    end

    @testset "Iteration through DNA Kmer" begin
        @test iterate(DNAKmer("ACTG")) == (DNA_A, 2)
        @test iterate(DNAKmer(""))     === nothing

        @test iterate(DNAKmer("ACTG"), 1) == (DNA_A, 2)
        @test iterate(DNAKmer("ACTG"), 4) == (DNA_G, 5)

        @test iterate(DNAKmer(""), 1)      === nothing
        @test iterate(DNAKmer("ACTG"), 1)  !== nothing
        @test iterate(DNAKmer("ACTG"), 4)  !== nothing
        @test iterate(DNAKmer("ACTG"), 5)  === nothing
        @test iterate(DNAKmer("ACTG"), -1) !== nothing

        dna_vec = [DNA_A, DNA_C, DNA_T, DNA_G]
        @test all([nt === dna_vec[i] for (i, nt) in enumerate(dna_kmer)])
    end

    @testset "Access RNA Kmer" begin
        @test rna_kmer[1] == RNA_A
        @test rna_kmer[2] == RNA_C
        @test rna_kmer[3] == RNA_U
        @test rna_kmer[4] == RNA_G

        # Access indexes out of bounds
        @test_throws Exception rna_kmer[-1]
        @test_throws Exception rna_kmer[0]
        @test_throws Exception rna_kmer[5]
        @test_throws Exception getindex(rna_kmer, -1)
        @test_throws Exception getindex(rna_kmer, 0)
        @test_throws Exception getindex(rna_kmer, 5)
    end

    @testset "Iteration through RNA Kmer" begin
        @test iterate(RNAKmer("ACUG")) == (RNA_A, 2)
        @test iterate(RNAKmer("")) === nothing

        @test iterate(RNAKmer("ACUG"), 1) == (RNA_A, 2)
        @test iterate(RNAKmer("ACUG"), 4) == (RNA_G, 5)

        @test iterate(RNAKmer(""), 1)      === nothing
        @test iterate(RNAKmer("ACUG"), 1)  !== nothing
        @test iterate(RNAKmer("ACUG"), 4)  !== nothing
        @test iterate(RNAKmer("ACUG"), 5)  === nothing
        @test iterate(RNAKmer("ACUG"), -1) !== nothing

        rna_vec = [RNA_A, RNA_C, RNA_U, RNA_G]
        @test all([nt === rna_vec[i] for (i, nt) in enumerate(rna_kmer)])
    end
end


text:
using Test, CommonLH, StructLH, CollegeStratCollege

function tf_test(switches)
    tf = typeof(switches);
    @testset "Tuition function $tf" begin
        StructLH.describe(switches)
        tf = init_tuition_fct(ObjectId(:test), switches);
        @test validate_tf(tf)

        t1 = tuition(tf, 2, 0.5, 0.6, 2);
        @test check_float(t1)
        t2 = tuition(tf, 3, [0.5, 0.6], 0.1, 3);
        @test check_float_array(t2, -1e5, 1e5)
        @test size(t2) == (2,)

        iCollege = CollInt(1);
        t3 = tuition(tf, iCollege, 0.5, 0.6, 3);
        if CollegeStratCollege.isfree(tf, iCollege)
            @test t3 == 0.0
        else
            @test t3 != 0.0
        end

        iCollege = CollInt(2);
        make_free!(switches, iCollege);
        tf = init_tuition_fct(ObjectId(:test), switches);
        t4 = tuition(tf, iCollege, 0.5, 0.6, 3);
        @test t4 == 0.0
	end
end

@testset "Tuition function" begin
    nc = 4;
    for switches in (test_tuition_by_qual_switches(nc),
        test_tuition_linear_switches(nc),
        test_tuition_linear_switches(nc, freeIdxV = [CollInt(1)]))

        tf_test(switches);
    end
end

# -------------


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Statistics, Test

tet4_meshfile = "test_problems_contact_3d/tet4.inp"
tet10_meshfile = "test_problems_contact_3d/tet10.inp"

function get_model(meshfile)
    mesh = abaqus_read_mesh(meshfile)

    upper = Problem(Elasticity, "UPPER", 3)
    upper_elements = create_elements(mesh, "UPPER")
    update!(upper_elements, "youngs modulus", 3*288.0)
    update!(upper_elements, "poissons ratio", 1/3)
    add_elements!(upper, upper_elements)

    lower = Problem(Elasticity, "LOWER", 3)
    lower_elements = create_elements(mesh, "LOWER")
    update!(lower_elements, "youngs modulus", 288.0)
    update!(lower_elements, "poissons ratio", 1/3)
    add_elements!(lower, lower_elements)

    bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
    bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
    update!(bc_upper_elements, "displacement 3", -0.4)
    add_elements!(bc_upper, bc_upper_elements)

    bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
    bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
    update!(bc_lower_elements, "displacement 3", 0.0)
    add_elements!(bc_lower, bc_lower_elements)

    # point-wise boundary conditions to prevent free body move
    nid1 = find_nearest_nodes(mesh, [0.0, 0.0, 0.0])[1]
    nid2 = find_nearest_nodes(mesh, [1.0, 0.0, 0.0])[1]
    nid3 = find_nearest_nodes(mesh, [0.0, 1.0, 0.0])[1]
    nid4 = find_nearest_nodes(mesh, [0.0, 0.0, 1.0])[1]
    nid5 = find_nearest_nodes(mesh, [1.0, 0.0, 1.0])[1]
    nid6 = find_nearest_nodes(mesh, [0.0, 1.0, 1.0])[1]

    bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
    # nodes in X2=0 plane
    bc_sym13_elements = [Element(Poi1, [j]) for j in [nid1, nid2, nid4, nid5]]
    update!(bc_sym13_elements, "geometry", mesh.nodes)
    update!(bc_sym13_elements, "displacement 2", 0.0)
    add_elements!(bc_sym13, bc_sym13_elements)

    bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
    # nodes in X1=0 plane
    bc_sym23_elements = [Element(Poi1, [j]) for j in [nid1, nid3, nid4, nid6]]
    update!(bc_sym23_elements, "geometry", mesh.nodes)
    update!(bc_sym23_elements, "displacement 1", 0.0)
    add_elements!(bc_sym23, bc_sym23_elements)

    interface = Problem(Contact, "LOWER_TO_UPPER", 3, "displacement")
    interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
    interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
    update!(interface_slave_elements, "master elements", interface_master_elements)
    add_elements!(interface, interface_slave_elements, interface_master_elements)
    interface.properties.contact_state_in_first_iteration = :AUTO

    #append!(bc_sym13.assembly.removed_dofs, [1316, 1319, 1388, 1358])
    #append!(bc_sym23.assembly.removed_dofs, [1492, 1627, 1387, 1597])
    #append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])

    analysis = Analysis(Nonlinear)
    add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
#    solver.properties.max_iterations = 5

    return analysis
end

## small sliding contact patch test, tet4 + standard basis
analysis = get_model(tet4_meshfile)
xdmf = Xdmf("contact_sl_lin_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
interface = get_problem(analysis, "LOWER_TO_UPPER")
interface.properties.dual_basis = false
run!(analysis)
u = interface("displacement", 0.0)
X = interface("geometry", 0.0)
# test postprocess of fields
postprocess!(interface, 0.0, Val{Symbol("contact pressure")})
contact_pressure = interface("contact pressure", 0.0)
@debug("Contact pressure in interface", contact_pressure)
u3 = [u[3] for u in values(u)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
maxpres = maximum(values(contact_pressure))
stdpres = std(values(contact_pressure))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@debug("max(contact_pressure) = $maxpres, std(contact_pressure) = $stdpres")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)
@test isapprox(maxpres, 172.8; atol=1.0e-6)

## small sliding contact patch test, tet4 + dual basis
analysis = get_model(tet4_meshfile)
xdmf = Xdmf("contact_dl_lin_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
interface = get_problem(analysis, "LOWER_TO_UPPER")
interface.properties.dual_basis = true
run!(analysis)
u = interface("displacement", 0.0)
X = interface("geometry", 0.0)
u3 = [u[3] for u in values(u)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)

## small sliding contact patch test, tet10 + standard basis
analysis = get_model(tet10_meshfile)
xdmf = Xdmf("contact_sl_quad_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
interface = get_problem(analysis, "LOWER_TO_UPPER")
interface.properties.dual_basis = false
run!(analysis)
u = interface("displacement", 0.0)
X = interface("geometry", 0.0)
u3 = [u[3] for u in values(u)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)

## small sliding contact patch test, tet10 + dual basis, alpha=0.2
analysis = get_model(tet10_meshfile)
xdmf = Xdmf("contact_dl_quad_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
interface = get_problem(analysis, "LOWER_TO_UPPER")
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
run!(analysis)
u = interface("displacement", 0.0)
X = interface("geometry", 0.0)
u3 = [u[3] for u in values(u)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)


text:
#############################################################################
##
##  SingerAlg.jl
##
##  Examine the Loewy structure of the algebras ``A[q,z] = A(q,n,e)``.
##

"""
  SingerAlg.jl provides Julia code for computations concerning the algebras
  ``A[q,z]`` that are studied in the following papers.

  [BHHK20a](https://doi.org/10.1016/j.jalgebra.2019.05.004)
  Breuer, T., Héthelyi, L., Horváth, E. and Külshammer, B.,
  The Loewy Structure of Certain Fixpoint Algebras, Part I,
  J. Algebra 558 (2020), 199–220.

  [BHHK20b](http://arxiv.org/abs/1912.03065v1)
  Breuer, T., Héthelyi, L., Horváth, E. and Külshammer, B.,
  The Loewy Structure of Certain Fixpoint Algebras, Part II,
  http://arxiv.org/abs/1912.03065v1.

  The Julia functions are used by the
  [GAP](https://www.gap-system.org) package
  [SingerAlg](http://www.math.rwth-aachen.de/~Thomas.Breuer/singeralg).
"""
module SingerAlg

##  One reason for importing `GAP.jl` is that we want to set some
##  GAP attributes from the Julia code as soon as the values become known.
##  Another reason is that we want to use GAP's `Factors` function
##  for decomposing integers into primes.
import GAP

## The factorization functions in Julia's `Primes` package (at least up to
## version 0.4.0) are substantially slower than the ones in GAP.
## We use the relevant GAP functions via GAP.jl.
include( "PrimesGAP.jl" )

include( "utils.jl" )
include( "loewy.jl" )
include( "mindeg.jl" )
include( "invar.jl" )

end

max_stars_repo_name: oscar-system/SingerAlg

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: oscar-system/SingerAlg

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: oscar-system/SingerAlg

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
module ProjectTemplate

using DifferentialEquations
using RecipesBase
using DiffEqBiological
using BlackBoxOptim
using CostFunctions
using DiffEqParameters

include("data.jl")
include("models.jl")
include("simulate.jl")
include("optimise.jl")
include("cost.jl")

Model = Union{DiffEqBase.AbstractReactionNetwork}

@param Param cost::Float64 metadata::Dict{Symbol, Any}
model(p::Param) = p.model
cost(p::Param) = p.cost

export simulate, optimise, cost, Param

end


text:
import Plots._current_plots_version

# Taken from MakieGallery
"""
Downloads the reference images from ReferenceImages for a specific version
"""
function download_reference(version = v"0.0.1")
    download_dir = abspath(@__DIR__, "reference_images")
    isdir(download_dir) || mkpath(download_dir)
    tarfile = joinpath(download_dir, "reference_images.zip")
    url = "https://github.com/JuliaPlots/PlotReferenceImages.jl/archive/v$(version).tar.gz"
    refpath = joinpath(download_dir, "PlotReferenceImages.jl-$(version)")
    if !isdir(refpath) # if not yet downloaded
        @info "downloading reference images for version $version"
        download(url, tarfile)
        BinaryProvider.unpack(tarfile, download_dir)
        # check again after download
        if !isdir(refpath)
            error("Something went wrong while downloading reference images. Plots can't be compared to references")
        else
            rm(tarfile, force = true)
        end
    else
        @info "using reference images for version $version (already downloaded)"
    end
    refpath
end

const ref_image_dir = download_reference()

function image_comparison_tests(pkg::Symbol, idx::Int; debug = false, popup = isinteractive(), sigma = [1,1], tol = 1e-2)
    Plots._debugMode.on = debug
    example = Plots._examples[idx]
    Plots.theme(:default)
    @info("Testing plot: $pkg:$idx:$(example.header)")
    backend(pkg)
    backend()
    default(size=(500,300))
    # ensure consistent results
    Random.seed!(1234)

    # reference image directory setup
    refdir = joinpath(ref_image_dir, "Plots", string(pkg))
    fn = "ref$idx.png"

    # firgure out version info
    vns = filter(x->x[1] != '.', readdir(refdir))
    versions = sort(VersionNumber.(vns), rev = true)
    versions = filter(v -> v <= _current_plots_version, versions)
    # @show refdir fn versions

    newdir = joinpath(refdir, string(_current_plots_version))
    newfn = joinpath(newdir, fn)

    # figure out which reference file we should compare to, by finding the highest versioned file
    reffn = nothing
    for v in versions
        tmpfn = joinpath(refdir, string(v), fn)
        if isfile(tmpfn)
            reffn = tmpfn
            break
        end
    end

    # now we have the fn (if any)... do the comparison
    # @show reffn
    if reffn == nothing
        reffn = newfn
    end
    # @show reffn
    # return

    # test function
    func = (fn, idx) -> begin
        expr = Expr(:block)
        append!(expr.args, example.exprs)
        eval(expr)
        png(fn)
    end

    # try
    #     run(`mkdir -p $newdir`)
    # catch err
    #     display(err)
    # end
    # # reffn = joinpath(refdir, "ref$idx.png")

    # the test
    vtest = VisualTest(func, reffn, idx)
    test_images(vtest, popup=popup, sigma=sigma, tol=tol, newfn = newfn)
end

function image_comparison_facts(pkg::Symbol;
                                skip = [],      # skip these examples (int index)
                                only = nothing, # limit to these examples (int index)
                                debug = false,  # print debug information?
                                sigma = [1,1],  # number of pixels to "blur"
                                tol = 1e-2)     # acceptable error (percent)
  for i in 1:length(Plots._examples)
    i in skip && continue
    if only == nothing || i in only
      @test image_comparison_tests(pkg, i, debug=debug, sigma=sigma, tol=tol) |> success == true
    end
  end
end


text:
using MomentTensors, Test

@testset "Operations" begin
    let m1 = MT(11, 22, 33, 12, 13, 23), m2 = MT(-11, -22, -33, -12, -13, -23)
        @test m1 + m2 == MT(0, 0, 0, 0, 0, 0)
        @test m1 == -m2
        @test -m1 == m2
        @test m1 + 1 == MT(12, 23, 34, 13, 14, 24)
        @test m1 - 1 == MT(10, 21, 32, 11, 12, 22)
        @test m1 + 3 == 3 + m1
        @test m1 - 4 == -(4 - m1)
        @test 2m1 == MT(22, 44, 66, 24, 26, 46)
        @test 2m1 == m1*2
        @test 1/m1 == MT(1 ./ (11, 22, 33, 12, 13, 23)...)
        @test m1/2 == MT(5.5, 11, 16.5, 6, 6.5, 11.5)
    end
end


text:
#### CONSEQUENT ###########
"
Structure that keeps the information related to one of the consequents of a production rule.

# Fields
- `_symbols::Array{UInt16}`: list with the symbols (in order) that represents the consequent.
- `_depth::Int16`: minimum depth that the generate tree would have, considering depth as the maximum
    of the minimum depths of the symbols.
- `_minProductions::Int16`: the number of productions that must be mandatory done through this consequent.
- `_recursive::Bool`: whether the consequent is recursive or not, meaning that the consequent
    has a symbol that is equal to its antecedent or there is production rule derivated that
    fulfill this statement.
- `_nSymbols::UInt8`: number of symbols in the consequent.

See also: [`Grammar`](@ref), [`ProductionRule`](@ref)
"
mutable struct Consequent
    _symbols::Array{UInt16}
    _depth::Int16
    _minProductions::Int16
    _recursive::Bool
    _nSymbols::UInt8

    Consequent(symbols::Array{UInt16}) = new(symbols, -32000, 32000, false, length(symbols))
    Consequent(symbols::Array{UInt16}, depth::Int16) = new(symbols, depth, 32000, false, length(symbols))
end

"""
    getNSymbols(genotype::Genotype)::UInt8

Returns the number of symbol of the consequent.
"""
getNSymbols(consequent::Consequent)::UInt8 = consequent._nSymbols

"""
    getDepth(consequent::Consequent)::Int16

Returns the the maximum of the minimum depths of the symbols.
"""
getDepth(consequent::Consequent)::Int16 = consequent._depth

"""
    isRecursive(consequent::Consequent)::Bool

Returns whether is recursive or not.
"""
isRecursive(consequent::Consequent)::Bool = consequent._recursive

"""
    getMinProductions(consequent::Consequent)::Int16

Returns the number of productions that must be mandatory done through this consequent.
"""
getMinProductions(consequent::Consequent)::Int16 = consequent._minProductions

#### PRODUCTION RULE ###########

"
Structure that keeps the information related to a production rule that is identified by antecedent.

# Fields
- `_antecedent::UInt16`: Non-terminal symbol.
- `_consequent::Array{Consequent}`
- `_depth::Int16`: depth that represent the minimum of consequents depth.
- `_nConsequent::UInt8`: number of consequents.

See also: [`Grammar`](@ref), [`Consequent`](@ref)
"
mutable struct ProductionRule
    _antecedent::UInt16
    _consequent::Array{Consequent}
    _depth::Int16
    _nConsequent::UInt8

    ProductionRule(antecedent::Integer, consequent::Array{Consequent}) = new(antecedent, consequent, -32000, length(consequent))
end

"""
    hasDepth(productionRule::ProductionRule)::Bool

Returns whether the depth is already set or not.
"""
hasDepth(productionRule::ProductionRule)::Bool = productionRule._depth != -32000

"""
    getNConsequent(productionRule::ProductionRule)::UInt8

Returns the number of consequent that the production rule has.
"""
getNConsequent(productionRule::ProductionRule)::UInt8 = productionRule._nConsequent

"""
    getNRules(productionRule::ProductionRule)::UInt8

Returns the number of consequent that the production rule has.
"""
getNRules(productionRule::ProductionRule)::UInt8 = getNConsequent(productionRule)

"""
    getDepth(productionRule::ProductionRule)::Int16

Returns the depth that represent the minimum of consequents depth.
"""
getDepth(productionRule::ProductionRule)::Int16 = productionRule._depth

"""
    getConsequent(productionRule::ProductionRule)::Consequent

Returns the consequent associated with index.
"""
getConsequent(productionRule::ProductionRule, index::UInt8)::Consequent = productionRule._consequent[index]

"""
    getPosibleConsequentProductions(productionRule::ProductionRule, productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}

Obtains the consequents indexes that have less or equal minimum productions than the productions specified.
"""
function getPosibleConsequentProductions(productionRule::ProductionRule, productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}
    index = UInt8[]
    count::UInt8 = 0x0
    for i=eachindex(productionRule._consequent)
        consequent = productionRule._consequent[i]
        if getMinProductions(consequent) <= productions
            push!(index, i)
            count += 0x1
        end
    end
    return index, count
end

"""
    getPosibleConsequentProductions(productionRule::ProductionRule, productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}

Obtains the consequents indexes that have less or equal minimum productions than the productions specified and are recursives.
"""
function getPosibleRecursiveProductions(productionRule::ProductionRule, productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}
    index = UInt8[]
    count::UInt8 = 0x0
    for i=eachindex(productionRule._consequent)
        consequent = productionRule._consequent[i]
        if getMinProductions(consequent) <= productions && isRecursive(consequent)
            push!(index, i)
            count += 0x1
        end
    end

    if isempty(index)
        index, count = getPosibleConsequentProductions(productionRule, productions)
    end
    return index, count
end

"""
    getPosibleConsequentDepth(productionRule::ProductionRule, depth::UInt8=0xffff, productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}

Obtains the consequents indexes that have less or equal minimum productions than the productions specified and
have less minDepth than the depth specified.
"""
function getPosibleConsequentDepth(productionRule::ProductionRule, depth::UInt8=0xffff, productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}
    index = UInt8[]
    count::UInt8 = 0x0
    for i=eachindex(productionRule._consequent)
        consequent = productionRule._consequent[i]
        if getDepth(consequent) < depth && getMinProductions(consequent) <= productions
            push!(index, i)
            count += 0x1
        end
    end
    return index, count
end

"""
    getPosibleRecursiveDepth(productionRule::ProductionRule, depth::UInt8=0xffff, productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}

Obtains the consequents indexes that have less or equal minimum productions than the productions specified,
have less minDepth than the depth specified and are recursive.
"""
function getPosibleRecursiveDepth(productionRule::ProductionRule, depth::UInt8=0xff, productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}
    index = UInt8[]
    count::UInt8 = 0x0
    for i=eachindex(productionRule._consequent)
        consequent = productionRule._consequent[i]
        if (getDepth(consequent) < depth
            && isRecursive(consequent)
            && getMinProductions(consequent) <= productions)
            push!(index, i)
            count += 0x1
        end
    end
    if isempty(index)
        index, count = getPosibleConsequentDepth(productionRule, depth, productions)
    end
    return index, count
end

#### GRAMMAR #########
"
Structure that keeps the information related to the grammar.

# Fields
- `_N::Array{String}`: non-terminal set.
- `_T::Array{String}`: terminal set.
- `_P::Array{ProductionRule}`: storages all the info about the production rules.
- `_S::String`: initial symbol of the grammar.
- `_nProductionRules::UInt16`: number of production rules

!!! note
    **Symbols** are `UInt`, that represent an index of the terminal array if
    symbol > 1000 or an index of the non terminal if symbol <= 1000
See also: [`Consequent`](@ref), [`ProductionRule`](@ref)
"
struct Grammar
    _N::Array{String}
    _T::Array{String}
    _P::Array{ProductionRule}
    _S::UInt16
    _nProductionRules::UInt16

    Grammar(N, T, P, S) = new(N, T, P, convert(UInt16, findfirst(x->x==S, N)), length(P))
end # struct

"""
    getInitialSymbol(gramm::Grammar)::UInt16

Returns the first symbol of the grammar.
"""
getInitialSymbol(gramm::Grammar)::UInt16 = gramm._S

"""
    ggetProductionRule(gramm::Grammar, id::UInt16)::ProductionRule

Obtains the production rule assoacited to id, that in this context is the same as
non terminal symbol.
"""
getProductionRule(gramm::Grammar, id::UInt16)::ProductionRule = gramm._P[id]

"""
    getDepth(gramm::Grammar, id::UInt16)::Int16

A interface for `getDepth(productionRule::ProductionRule)`.
"""
getDepth(gramm::Grammar, id::UInt16)::Int16 = getDepth(getProductionRule(gramm, id))

"""
    getRepresentation(gramm::Grammar, id::UInt16)::String

Returns the correspondence of a numeric symbol, a string.
"""
getRepresentation(gramm::Grammar, id::UInt16)::String = isTerminalSymbol(gramm, id) ? gramm._T[id-0x03e8] : gramm._N[id]

"""
    isTerminalSymbol(gramm::Grammar, id::UInt16)::Bool

Returns the correspondence of a numeric symbol, a string.
"""
isTerminalSymbol(gramm::Grammar, id::UInt16)::Bool = id > 0x03e7

"""
    getNRules(gramm::Grammar)::UInt16

Returns the number of rules in the grammar, meaning all the consequent.
"""
function getNRules(gramm::Grammar)::UInt16
    nRules::UInt16 = 0x0
    for i=eachindex(gramm._P)
        nRules = nRules + getNRules(gramm._P[i])
    end
    return nRules
end

"""
    getNProductionRules(gramm::Grammar)::UInt16

Returns the number of production rules in the grammar.
"""
getNProductionRules(gramm::Grammar)::UInt16 = gramm._nProductionRules

"""
    getNConsequent(gramm::Grammar, id::UInt16)::UInt8

A interface for `getNConsequent(productionRule::ProductionRule)`.
"""
getNConsequent(gramm::Grammar, id::UInt16)::UInt8 = getNConsequent(getProductionRule(gramm, id))

"""
    getPosibleConsequentProductions(gramm::Grammar, antecedent::UInt16,
            productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}

A interface for `getPosibleConsequentProductions(productionRule::ProductionRule, productions)`.
"""
getPosibleConsequentProductions(gramm::Grammar, antecedent::UInt16,
            productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8} = getPosibleConsequentProductions(getProductionRule(gramm, antecedent), productions)

"""
    getPosibleRecursiveProductions(gramm::Grammar, antecedent::UInt16,
            productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}

A interface for `getPosibleRecursiveProductions(productionRule::ProductionRule, productions)`.
"""
getPosibleRecursiveProductions(gramm::Grammar, antecedent::UInt16,
        productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8} = getPosibleRecursiveProductions(getProductionRule(gramm, antecedent), productions)

"""
    getPosibleConsequentDepth(gramm::Grammar, antecedent::UInt16, depth::UInt8=0xff,
            productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}

A interface for `getPosibleConsequentDepth(productionRule::ProductionRule, depth, productions)`.
"""
getPosibleConsequentDepth(gramm::Grammar, antecedent::UInt16, depth::UInt8=0xff,
        productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8} = getPosibleConsequentDepth(getProductionRule(gramm, antecedent), depth, productions)

"""
    getPosibleRecursiveDepth(gramm::Grammar, antecedent::UInt16, depth::UInt8=0xff,
            productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8}

A interface for `getPosibleRecursiveDepth(productionRule::ProductionRule, depth, productions)`.
"""
getPosibleRecursiveDepth(gramm::Grammar, antecedent::UInt16, depth::UInt8=0xff,
        productions::UInt16=0xffff)::Tuple{Array{UInt8}, UInt8} = getPosibleRecursiveDepth(getProductionRule(gramm, antecedent), depth, productions)

"""
    getConsequent(gramm::Grammar, antecedent::UInt16, index::UInt8)::Consequent

A interface for `getConsequent(productionRule::ProductionRule, index)`.
"""
getConsequent(gramm::Grammar, antecedent::UInt16, index::UInt8)::Consequent = getConsequent(getProductionRule(gramm, antecedent), index)

"""
    getMaxSymbols(gramm::Grammar)

Returns the number of non terminal symbols of the rule that has more non terminal symbol.
"""
function getMaxSymbols(gramm::Grammar)
    max = 0
    for productionRule in gramm._P
        for consequent in productionRule._consequent
            current = 0
            for i in eachindex(consequent._symbols)
                if !isTerminalSymbol(gramm, consequent._symbols[i])
                    current += 1
                end
            end
            if current > max
                max = current
            end
        end
    end
    return max
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# Nested distance optimization with Optim for large trees

using ScenTreesMakie
using Optim
using JLD

Tree0 = load("Trees/311000/binary.jld","Tree0");

#############################
# Regularization Ansatz

function myfunc(x)

    trr1.state = x # update states Vector
    #trr1.state[1] = trr2.state[1]
    dist = nestedWasserstein(trr1,trr2,2.0).^2.0

    return dist

end

###################
function myConstraint(x)

    # lower bound such that process is monotonically increasing
    lb = trr1.state[trr1.parent[2:end]]
    tmp = sum( max.([0;lb] .- x,0))

    #upper bound such that fitted tree does not exceed original tree (significantly)
    ub = maximum(trr2.state);
    tmp2 = sum(max.(x .- ub,0) )

    return( tmp + tmp2 )
end

function myRegularizedFunction(x)
    return( myfunc(x) + 1000*myConstraint(x) )

end

bla =[1,2,4,8,16,33,66,133,267,535,1070,2141,4282,8565,17131,34263,68527] # from other laptop; works only with CBDX_simple_binary
plot(Tree0)

# Save list of trees in stage
t = 6; trrs = [];
#for i=1:length(nodes(Tree0,stage))
#    println(i)
#    push!(trrs,part_tree(Tree0,nodes(Tree0,stage)[i]  ))
#end
#save("partialTrees06-CBDXsimpleBinary.jld","trrs",trrs)

initialTree = deepcopy(Tree0);
lastNode = nodes(Tree0,t-1)[end]
initialTree.state=initialTree.state[1:lastNode];
initialTree.probability=initialTree.probability[1:lastNode];
initialTree.parent=initialTree.parent[1:lastNode];
initialTree.children=initialTree.children[1:lastNode];

trrs= load("partialTrees06-CBDXsimpleBinary.jld","trrs");
for i=1:(length(trrs)-63)
    global trr1 = deepcopy(trrs[i]);
    global trr2 = deepcopy(trrs[i]);
    n = length(trr1.state);
    trr1.state[2:end] = max.(trr1.state[trr1.parent[2:end]] .- trr1.state[2:end],0) + trr1.state[2:end];
    x0 = trr1.state;

    result = Optim.optimize(myRegularizedFunction, x0, ParticleSwarm(n_particles = 5), Optim.Options(iterations=200, show_trace=true)  )
#    result = Optim.optimize(myRegularizedFunction, x0, NelderMead(), Optim.Options(g_tol = 1e-8, iterations=500, show_trace=true)  )
    trrs[i] = trr1;
end

target = deepcopy(Tree0)
target = combine_tree(target,initialTree,trrs)

result = Optim.optimize(myRegularizedFunction, x0, ParticleSwarm(n_particles = 5), Optim.Options(iterations=200, show_trace=true)  )

#result = Optim.optimize(myRegularizedFunction, x0, NelderMead(), Optim.Options(g_tol = 1e-8, iterations=500, show_trace=true)  )
Optim.minimizer(result)

plot(trr1,nothing,true)
plot(trr2,nothing,true)


text:
# See: Neutral models of land-use change as benchmarks in the assessment of model performance
# Alex Hagen-Zanker and Gilles Lajoie

function category_list(A, cat)
    list = PreallocatedUnorderedList(CartesianIndex{2}; maxlength=length(A))
    for I in CartesianIndices(A)
        if A[I] == cat
            push!(list, I)
        end
    end
    return list
end

adjust(A, B, categories) = map(cat -> adjust(A, B, cat), categories)
adjust(A, B, cat::Int) = count(==(cat), A) - count(==(cat), B)

abstract type NeutralLUCModel end
struct RandomConstraintMatch <: NeutralLUCModel end
struct GrowingClusters  <: NeutralLUCModel end

function sim(::RandomConstraintMatch, init, target, categories)
    result = copy(init)
    catlists = map(c -> category_list(init, c), categories)
    celllist = PreallocatedUnorderedList(CartesianIndex{2}; maxlength=length(init))
    adj = adjust(init, target, categories)
    for cat in categories
        adj[cat] > 0 || continue
        catlist = catlists[cat]
        for sample in 1:adj[cat]
            read_random!(celllist, catlist, cat)
        end
    end
    for cat in categories
        adj[cat] < 0 || continue
        for sample in 1:-adj[cat]
            write_random!(result, celllist, cat)
        end
    end
    return result
end

function read_random!(celllist, catlist, cat)
    i = rand(1:length(catlist))
    push!(celllist, catlist[i])
    deleteat!(catlist, i)
end

function write_random!(result, celllist, cat)
    i = rand(1:length(celllist))
    result[celllist[i]] = cat
    deleteat!(celllist, i)
end

function sim(::GrowingClusters, init, target, categories; neighborhood=VonNeumann{1}())
    result = copy(init)
    ax = unpad_axes(init, neighborhood)
    inds = CartesianIndices(ax)
    edgetype = NamedTuple{(:cell, :cat, :hoodcat),Tuple{CartesianIndex{2},Int,Int}}
    edgelist = UnorderedList(edgetype, length(neighborhood) * length(inds))
    adj = [map(cat -> adjust(view(init, ax...), view(target, ax...), cat), categories)...]
    while !all(adj .== 0)
        for cell in inds
            hood = unsafe_updatewindow(neighborhood, result, cell)
            cat = result[cell]
            add_hood_edges!(edgelist, hood, adj, cell, cat)
        end
        if isempty(edgelist)
            for cat in categories
                if adj[cat] < 0
                    while adj[randomcat] <= 0
                        cell = rand(inds)
                        randomcat = result[cell]
                    end
                    result[cell] = cat
                    adj[randomcat] -= 1
                    adj[cat] += 1
                end
            end
        end
        while !isempty(edgelist)
            i = rand(1:length(edgelist))
            edge = edgelist[i]
            deleteat!(edgelist, i)
            write_if_matches!(result, adj, edge)
        end
    end
    return result
end

# These are mostly separated as function barriers for type stability
function write_if_matches!(result, adj, edge)
    stillover = adj[edge.cat] > 0
    stillunder = adj[edge.hoodcat] < 0
    stillsame = result[edge.cell] == edge.cat
    if stillover && stillunder && stillsame
        adj[edge.cat] -= 1
        adj[edge.hoodcat] += 1
        result[edge.cell] = edge.hoodcat
    end
end

function add_hood_edges!(edgelist, hood, adj, cell, cat)
    for hoodcat in hood
        if (adj[cat] > 0) && (adj[hoodcat] < 0)
            push!(edgelist, (; cell, cat, hoodcat))
        end
    end
end


text:
###Policies that determine the options in the model evaluation
immutable LatencyCalculation <: GeneralizedMetropolisHastings.AbstractPolicyTrait
    trait::Symbol
    LatencyCalculation(s::Symbol) = (@assert in(s,[:deterministic,:lognormal]) ; new(s))
end

immutable RefractoryCalculation <: GeneralizedMetropolisHastings.AbstractPolicyTrait
    trait::Symbol
    RefractoryCalculation(s::Symbol) = (@assert in(s,[:deterministic,:lognormal]) ; new(s))
end

immutable BumpShape <: GeneralizedMetropolisHastings.AbstractPolicyTrait
    trait::Symbol
    BumpShape(s::Symbol) = (@assert in(s,[:fixed,:sample]) ; new(s))
end

immutable PlotMode <: GeneralizedMetropolisHastings.AbstractPolicyTrait
    trait::Symbol
    PlotMode(s::Symbol) = (@assert in(s,[:repl,:notebook]) ; new(s))
end

immutable FixedBumpShapeParameters{T<:AbstractFloat}
    amplitude::T
    shape::T
    scale::T
    cutoff::T
end

_traitname(t::LatencyCalculation) = :latency
_traitname(t::RefractoryCalculation) = :refractory
_traitname(t::BumpShape) = :bump
_traitname(t::PlotMode) = :plotmode

_trait(::Type{Val{:latency}},t::Symbol) = LatencyCalculation(t)
_trait(::Type{Val{:refractory}},t::Symbol) = RefractoryCalculation(t)
_trait(::Type{Val{:bump}},t::Symbol) = BumpShape(t)
_trait(::Type{Val{:plotmode}},t::Symbol) = PlotMode(t)

paramkeys(t::GeneralizedMetropolisHastings.AbstractPolicyTrait) = tuple([Symbol(_traitname(t),x) for x in _paramkeys(t)]...)
parampriors(t::GeneralizedMetropolisHastings.AbstractPolicyTrait) = _parampriors(t)

function _paramkeys(t::Union{LatencyCalculation,RefractoryCalculation})
    if t.trait == :deterministic
        return (:value,)
    elseif t.trait == :lognormal
        return (:location,:scale)
    else
        return ()
    end
end

function _paramkeys(t::BumpShape)
    if t.trait == :sample
        return (:amplitude,:shape,:scale)
    else
        return ()
    end
end

function _parampriors(t::Union{LatencyCalculation,RefractoryCalculation})
    if t.trait == :deterministic
        return (:uniform,)
    elseif t.trait == :lognormal
        return (:uniform,:uniform)
    else
        return ()
    end
end

function _parampriors(t::BumpShape)
    if t.trait == :sample
        return (:uniform,:lognormal,:lognormal)
    else
        return ()
    end
end

fixedbumpshape(;amplitude =4.0f0,shape =3.0f0,scale =2.5f0,cutoff =0.025f0) = FixedBumpShapeParameters(amplitude,shape,scale,cutoff)
fixedbumpvalues(f::FixedBumpShapeParameters) = (f.amplitude,f.shape,f.scale,f.cutoff)

immutable PhotoReceptorPolicy <: GeneralizedMetropolisHastings.AbstractPolicy
    latency::LatencyCalculation
    refractory::RefractoryCalculation
    bump::BumpShape
    fixedbump::FixedBumpShapeParameters
    PhotoReceptorPolicy(l,r,b,d) = new(l,r,b,d)
end

function _policy(::Type{Val{:photoreceptor}};latency =:lognormal,refractory =:lognormal,bump =:fixed,fixedbump =fixedbumpshape())
    PhotoReceptorPolicy(trait(:latency,latency),trait(:refractory,refractory),trait(:bump,bump),fixedbump)
end

function show(io::IO,b::FixedBumpShapeParameters)
    println(io,"  FixedBumpShapeParameters with values:")
    println(io,"    amplitude = ",b.amplitude)
    println(io,"    shape = ",b.shape)
    println(io,"    scale = ",b.scale)
    println(io,"    cutoff = ",b.cutoff)
    nothing
end

function show(io::IO,p::PhotoReceptorPolicy)
    println(io,"PhotoReceptorPolicy with traits:")
    println(io,"  latency = ",traitvalue(p.latency))
    println(io,"  refractory = ",traitvalue(p.refractory))
    println(io,"  bump = ",traitvalue(p.bump))
    show(io,p.fixedbump)
    nothing
end


text:
## --------------------------------------------------
## no Doc trick
Base.Docs.getdoc(u::Sym) = Base.Docs.getdoc(getpy(u))
## --------------------------------------------------

is_(k::Symbol, x::Sym) = is_(k, getpy(x))
function is_(k::Symbol, x::Py)::Bool
    key = "is_$(string(k))"
    val = pygetattr(x, key, false)
    Bool(val)
end

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
# if the environment is not defined, use this default
const GENIE_ENV = "dev"

if ! haskey(ENV, "GENIE_ENV")
  ENV["GENIE_ENV"] = GENIE_ENV
end
if ! haskey(ENV, "NWORKERS") && in("s", ARGS)
  ENV["NWORKERS"] = 1
end


text:
module glmSolverjl

include("io.jl")
include("link.jl")
include("family.jl")
include("solver.jl")
include("tools.jl")
include("fit.jl")
include("sample.jl")
include("simulate.jl")

export matrixToBlock,
       write2DArray,
       writeBlockMatrix,
       read2DArray,
       readBlockMatrix,
       writeNDArray,
       readNDArray,
       AbstractLink,
       IdentityLink,
       LogLink,
       InverseLink,
       NegativeBinomialLink,
       LogitLink,
       ProbitLink,
       CauchitLink,
       OddsPowerLink,
       LogComplimentLink,
       LogLogLink,
       ComplementaryLogLogLink,
       PowerLink,
       linkfun,
       deta_dmu,
       linkinv,
       AbstractDistribution,
       BernoulliDistribution,
       BinomialDistribution,
       GammaDistribution,
       PoissonDistribution,
       GaussianDistribution,
       InverseGaussianDistribution,
       NegativeBernoulliDistribution,
       PowerDistribution,
       init!,
       variance,
       devianceResiduals,
       AbstractSolver,
       VanillaSolver,
       QRSolver,
       solve,
       AbstractInverse,
       GETRIInverse,
       POTRIInverse,
       SYTRFInverse,
       GESVDInverse,
       inv,
       GESVSolver,
       POSVSolver,
       SYSVSolver,
       GELSSolver,
       GELSYSolver,
       GELSDSolver,
       AbstractGradientDescentSolver,
       GradientDescentSolver,
       MomentumSolver,
       NesterovSolver,
       AdagradSolver,
       AdadeltaSolver,
       RMSpropSolver,
       AdamSolver,
       AdaMaxSolver,
       NAdamSolver,
       AMSGradSolver,
       XWX,
       solve!,
       cov,
       Z,
       W,
       Control,
       absoluteError,
       relativeError,
       AbstractMatrixType,
       RegularData,
       AbstractGLM,
       GLM,
       GLMBlock1D,
       Block1D,
       Block1DParallel,
       glm,
       generateRandomMatrix,
       cov2cor,
       cor2cov,
       mvrnorm,
       AbstractSampleDistribution,
       BetaSampleDistribution,
       sample,
       UniformSampleDistribution,
       min,
       max,
       range,
       I,
       AbstractRandomCorrelationMatrix,
       BetaGenerator,
       OnionGenerator,
       UniformGenerator,
       VineGenerator,
       randomCorrelationMatrix,
       simulateData,
       AbstractPoissonDistribution,
       PoissonSampleDistribution

end # module

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
"""
Imputes gradients based on a vector of x and y coordinates.
### Takes
 * x - A Float64 vector of x coordinates
 * y - A Float64 vector of y coordinates

### Returns
 * A Float64 vector of gradients for each input point
"""
function imputeGradients(x,y)
     n = length(x)
     # Judd (1998), page 233, second last equation
     L = sqrt( (x[2:n]-x[1:(n-1)]).^2 + (y[2:n]-y[1:(n-1)]).^2)
     # Judd (1998), page 233, last equation
     d = (y[2:n]-y[1:(n-1)])./(x[2:n]-x[1:(n-1)])
     # Judd (1998), page 234, Eqn 6.11.6
     Conditionsi = d[1:(n-2)].*d[2:(n-1)] .> 0
     MiddleSiwithoutApplyingCondition = (L[1:(n-2)].*d[1:(n-2)]+L[2:(n-1)].* d[2:(n-1)]) ./ (L[1:(n-2)]+L[2:(n-1)])
     sb = Conditionsi .* MiddleSiwithoutApplyingCondition
     # Judd (1998), page 234, Second Equation line plus 6.11.6 gives this array of slopes.
     ff = [((-sb[1]+3*d[1])/2);  sb ;  ((3*d[n-1]-sb[n-2])/2)]
     return ff
 end

 """
 Creates a spline defined by interval starts IntStarts and quadratic coefficients SpCoefs which evaluates an input point.
### Takes
  * IntStarts - A Float64 vector that gives the starting points of intervals (in the x plane)
  * SpCoefs - A 3 column matrix with the same number of rows as the length of the IntStarts vector. The first column is the coefficient of the quadratic term, the second column for the linear term. The third column is the constant.

### Returns
  * A spline function that takes a single Float64 input and returns the spline value at that point.
 """
function ppmak(IntStarts,SpCoefs)
  function sp(PointToExamine)
    IntervalNum = searchsortedlast(IntStarts, PointToExamine)
    IntervalNum = max(IntervalNum, 1)
    xmt = PointToExamine - IntStarts[IntervalNum]
    Coefs = SpCoefs[ IntervalNum , :]
    return reshape(Coefs' * [xmt^2 xmt 1]', 1)[1]
  end
  function Vsp(PointToExamine)
    return map(x -> sp(x), PointToExamine)
  end
  return Vsp
end

 """
 Creates the derivative function of the spline defined by interval starts IntStarts and quadratic coefficients SpCoefs which evaluates an input point.
### Takes
  * IntStarts - A Float64 vector that gives the starting points of intervals (in the x plane)
  * SpCoefs - A 3 column matrix with the same number of rows as the length of the IntStarts vector. The first column is the coefficient of the quadratic term, the second column for the linear term. The third column is the constant.

### Returns
 * The derivative function that takes a single Float64 input and returns the derivative at that point.
 """
function ppmakDeriv(IntStarts,SpCoefs)
  function sp(PointToExamine)
    IntervalNum = searchsortedlast(IntStarts, PointToExamine)
    IntervalNum = max(IntervalNum, 1)
    xmt = PointToExamine - IntStarts[IntervalNum]
    Coefs = SpCoefs[ IntervalNum , :]
    return reshape(Coefs' * [2*xmt 1 0]', 1)[1]
  end
  function Vsp(PointToExamine)
    return map(x -> sp(x), PointToExamine)
  end
  return Vsp
end

"""
Creates the second derivative function of the spline defined by interval starts IntStarts and quadratic coefficients SpCoefs which evaluates an input point.
### Takes
 * IntStarts - A Float64 vector that gives the starting points of intervals (in the x plane)
 * SpCoefs - A 3 column matrix with the same number of rows as the length of the IntStarts vector. The first column is the coefficient of the quadratic term, the second column for the linear term. The third column is the constant.

### Returns
 * The second derivative function that takes a single Float64 input and returns the second derivative at that point.
"""
function ppmak2Deriv(IntStarts,SpCoefs)
  function sp(PointToExamine)
    IntervalNum = searchsortedlast(IntStarts, PointToExamine)
    IntervalNum = max(IntervalNum, 1)
    xmt = PointToExamine - IntStarts[IntervalNum]
    Coefs = SpCoefs[ IntervalNum , :]
    return reshape(Coefs' * [2 0 0]', 1)[1]
  end
  function Vsp(PointToExamine)
    return map(x -> sp(x), PointToExamine)
  end
  return Vsp
end

"""
Splits an interval into 2 subintervals and creates the quadratic coefficients
### Takes
 * s - A 2 entry Float64 vector with gradients at either end of the interval
 * z - A 2 entry Float64 vector with y values at either end of the interval
 * Smallt - A 2 entry Float64 vector with x values at either end of the interval

### Returns
 * A 2 x 5 matrix. The first column is the x values of start of the two subintervals. The second column is the ends. The last 3 columns are quadratic coefficients in two subintervals.
"""
function schumakerIndInterval(s,z,Smallt)
   # The SchumakerIndInterval function takes in each interval individually
   # and returns the location of the knot as well as the quadratic coefficients in each subinterval.

   # Judd (1998), page 232, Lemma 6.11.1 provides this if condition:
   if (sum(s)*(Smallt[2]-Smallt[1]) == 2*(z[2]-z[1]))
     tsi = Smallt[2]
   else
     # Judd (1998), page 233, Algorithm 6.3 along with equations 6.11.4 and 6.11.5 provide this whole section
     delta = (z[2] -z[1])/(Smallt[2]-Smallt[1])
     Condition = ((s[1]-delta)*(s[2]-delta) >= 0)
     Condition2 = abs(s[2]-delta) < abs(s[1]-delta)
     if (Condition)
       tsi = sum(Smallt)/2
     elseif (Condition2)
       tsi = (Smallt[1] + (Smallt[2]-Smallt[1])*(s[2]-delta)/(s[2]-s[1]))
     else
       tsi = (Smallt[2] + (Smallt[2]-Smallt[1])*(s[1]-delta)/(s[2]-s[1]))
     end
   end

   # Judd (1998), page 232, 3rd last equation of page.
   alpha = tsi-Smallt[1]
   beta = Smallt[2]-tsi
   # Judd (1998), page 232, 4th last equation of page.
   sbar = (2*(z[2]-z[1])-(alpha*s[1]+beta*s[2]))/(Smallt[2]-Smallt[1])
   # Judd (1998), page 232, 3rd equation of page. (C1, B1, A1)
   Coeffs1 = [ (sbar-s[1])/(2*alpha)  s[1]  z[1] ]
   if (beta == 0)
     Coeffs2 = Coeffs1
   else
     # Judd (1998), page 232, 4th equation of page. (C2, B2, A2)
     Coeffs2 = [ (s[2]-sbar)/(2*beta)  sbar  Coeffs1 * [alpha^2, alpha, 1] ]
   end
   Machine4Epsilon = 4*eps()
     if (tsi  <  Smallt[1] + Machine4Epsilon )
         return [Smallt[1] Smallt[2] Coeffs2]
     elseif (tsi + Machine4Epsilon > Smallt[2] )
         return [Smallt[1] Smallt[2] Coeffs1]
     else
         return [Smallt[1] tsi Coeffs1 ; tsi Smallt[2] Coeffs2]
     end
 end

 """
 Calls SchumakerIndInterval many times to get full set of spline intervals and coefficients. Then calls extrapolation for out of sample behaviour
### Takes
 * gradients - A Float64 vector of gradients at each point
 * x - A Float64 vector of x coordinates
 * y - A Float64 vector of y coordinates
 * extrapolation - A string in ("Curve", "Linear", "Constant") that gives behaviour outside of interpolation range.

### Returns
 * A vector of interval starts
 * A vector of interval ends
 * A matrix of all coefficients
  """
 function getCoefficientMatrix(gradients,y,x, extrapolation)
   n = length(x)
   fullMatrix = schumakerIndInterval([gradients[1] gradients[2]], [y[1] y[2]], [x[1] x[2]] )
    for intrval = 2:(n-1)
      Smallt = [ x[intrval] , x[intrval + 1] ]
      s = [ y[intrval], y[intrval + 1] ]
      z = [ gradients[intrval], gradients[intrval + 1] ]
      intMatrix = schumakerIndInterval(z,s,Smallt)
      fullMatrix = vcat(fullMatrix,intMatrix)
    end
    fullMatrix = extrapolate(fullMatrix, extrapolation, x, y)
   return fullMatrix[:,1], fullMatrix[:,2], fullMatrix[:,3:5]
 end

"""
 Adds a row on top and bottom of coefficient matrix to give out of sample prediction.
### Takes
 * fullMatrix - output from GetCoefficientMatrix first few lines
 * extrapolation - A string in ("Curve", "Linear", "Constant") that gives behaviour outside of interpolation range.
 * x - A Float64 vector of x coordinates
 * y - A Float64 vector of y coordinates

### Returns
  * A new version of fullMatrix with out of sample prediction built into it.
"""
function extrapolate(fullMatrix, extrapolation, x, y)

  if (extrapolation == "Curve")
    return fullMatrix
  end

  dim = size(fullMatrix)[1]

  Botx   = fullMatrix[1,1]
  Boty   = y[1]

  if (extrapolation == "Linear")
    BotB = fullMatrix[1 , 4]
    BotC   = Boty - BotB
  else
    BotB = 0.0
    BotC = Boty
  end

  BotRow = [ Botx-1, Botx, 0.0, BotB, BotC]

  Topx = fullMatrix[dim,2]
  Topy = y[length(y)]

  if (extrapolation == "Linear")
    TopB = fullMatrix[dim ,4]
    TopC = Topy
  else
    TopB = 0.0
    TopC = Topy
  end

  TopRow = [ Topx, Topx + 1, 0.0 ,TopB ,TopC]

  fullMatrix = vcat(BotRow' , fullMatrix,  TopRow')

  return fullMatrix
end

"""
Creates splines for a given set of x and y values (and optionally gradients) and the first and second derivatives of this spline.
### Takes
* x - A Float64 vector of x coordinates.
* y - A Float64 vector of y coordinates.
* gradients (optional)- A Float64 vector of gradients at each point. If not supplied these are imputed from x and y.
* extrapolation (optional) - This should be a string in ("Curve", "Linear", "Constant") specifying how to interpolate outside of the sample domain. By default it is "curve" which extends out the first and last quadratic curves. The other options are "Linear" which extends the line (from first and last curve) out from the first and last point and "Constant" which extends out the y value at the first and last point.

### Returns
* A spline which takes and input value and returns the spline y value.
* The derivative of this spline.
* The second derivative of this spline.
 """
function schumaker(x,y,gradients = "Not-Supplied", extrapolation = ("Curve", "Linear", "Constant"))
  # This is the main function of the package that creates and returns a schumaker spline.
  # Inputs : An x array, the corresponding y array and optionally the gradients at each point.
  #           If gradients are not input then they are estimated.
  # Outputs: A schumaker spline function. Its derivative and its second derivative.

  if (extrapolation == ("Curve", "Linear", "Constant"))
    extrapolation = "Curve"
  end

  if (gradients == "Not-Supplied")
     gradients = imputeGradients(x,y)
  end

  IntStarts, IntEnds, SpCoefs = getCoefficientMatrix(gradients,y,x, extrapolation)

  Sp      = ppmak(IntStarts,SpCoefs)
  SpDeriv = ppmakDeriv(IntStarts,SpCoefs)
  SpDeriv2= ppmak2Deriv(IntStarts,SpCoefs)
  return Sp, SpDeriv, SpDeriv2
end


text:
# Example: Compare the performance of some math functions

using BenchmarkLite

## Define VecMath procedures
#
# The configuration is the vector length, represented by an integer
#

type VecMath{Op} <: Proc end

type Sqrt end
calc(::Sqrt, x) = sqrt(x)

type Exp end
calc(::Exp, x) = exp(x)

type Log end
calc(::Log, x) = log(x)

type Sin end
calc(::Sin, x) = sin(x)

# procedure name
Base.string{Op}(::VecMath{Op}) = string("vec-", lowercase("$Op"))

# pre-allocated arrays for running the procedure
typealias FVecPair (Vector{Float64},Vector{Float64})

# procedure codes
Base.length(p::VecMath, n::Int) = n

Base.isvalid(p::VecMath, n::Int) = (n > 0)

Base.start(p::VecMath, n::Int) = (rand(n), zeros(n))

function Base.run{Op}(p::VecMath{Op}, n::Int, s::FVecPair)
    x, y = s
    op = Op()
    for i = 1:n
        @inbounds y[i] = calc(op, x[i])
    end
end

Base.done(p::VecMath, n, s) = nothing

## perform the benchmark

procs = Proc[ VecMath{Sqrt}(),
              VecMath{Exp}(),
              VecMath{Log}(),
              VecMath{Sin}() ]

cfgs = 2 .^ (4:12)

println("Running log:")
println("--------------------")
rtable = run(procs, cfgs; duration=0.2)
@assert isa(rtable, BenchmarkTable)
println()

## show results

show(rtable; unit=:mps, cfghead="len")

# open("mathfuns.bench.csv", "w") do f
#   writecsv(f, rtable)
# end


text:
module ContactAlgorithms

include("constraints.jl")
export c_wall, c_wall_friction

end

max_stars_repo_name: SteffenPL/ContactAlgorithms

max_issues_repo_name: SteffenPL/ContactAlgorithms


text:
abstract type DescentDirectionAllocs{T <: Real, S} end

mutable struct DescentDirectionAllocsPC{T <: Real, S} <: DescentDirectionAllocs{T, S}
  Δxy_aff::S # affine-step solution of the augmented system [Δx_aff; Δy_aff], size nvar + ncon
  Δs_l_aff::S # size nlow
  Δs_u_aff::S # size nupp
  x_m_l_αΔ_aff::S # x + α_aff * Δxy_aff - lvar , size nlow
  u_m_x_αΔ_aff::S # uvar - (x + α_aff * Δxy_aff) , size nupp
  s_l_αΔ_aff::S # s_l + α_aff * Δs_l_aff , size nlow
  s_u_αΔ_aff::S # s_u + α_aff * Δs_u_aff , size nupp
  rxs_l::S # - σ * μ * e + ΔX_aff * Δ_S_l_aff , size nlow
  rxs_u::S # σ * μ * e + ΔX_aff * Δ_S_u_aff , size nupp
  function DescentDirectionAllocsPC(
    Δxy_aff::AbstractVector{T},
    Δs_l_aff::AbstractVector{T},
    Δs_u_aff::AbstractVector{T},
    x_m_l_αΔ_aff::AbstractVector{T},
    u_m_x_αΔ_aff::AbstractVector{T},
    s_l_αΔ_aff::AbstractVector{T},
    s_u_αΔ_aff::AbstractVector{T},
    rxs_l::AbstractVector{T},
    rxs_u::AbstractVector{T},
  ) where {T <: Real}
    S = typeof(Δxy_aff)
    return new{T, S}(
      Δxy_aff,
      Δs_l_aff,
      Δs_u_aff,
      x_m_l_αΔ_aff,
      u_m_x_αΔ_aff,
      s_l_αΔ_aff,
      s_u_αΔ_aff,
      rxs_l,
      rxs_u,
    )
  end
end

DescentDirectionAllocsPC(id::QM_IntData, S::DataType) where {T <: Real} = DescentDirectionAllocsPC(
  S(undef, id.nvar + id.ncon), # Δxy_aff
  S(undef, id.nlow), # Δs_l_aff
  S(undef, id.nupp), # Δs_u_aff
  S(undef, id.nlow), # x_m_l_αΔ_aff
  S(undef, id.nupp), # u_m_x_αΔ_aff
  S(undef, id.nlow), # s_l_αΔ_aff
  S(undef, id.nupp), # s_u_αΔ_aff
  S(undef, id.nlow), # rxs_l
  S(undef, id.nupp),  # rxs_u
)

convert(
  ::Type{<:DescentDirectionAllocs{T, S}},
  dda::DescentDirectionAllocsPC{T0, S0},
) where {T <: Real, S, T0 <: Real, S0} = DescentDirectionAllocsPC(
  convert(S.name.wrapper{T, 1}, dda.Δxy_aff),
  convert(S.name.wrapper{T, 1}, dda.Δs_l_aff),
  convert(S.name.wrapper{T, 1}, dda.Δs_u_aff),
  convert(S.name.wrapper{T, 1}, dda.x_m_l_αΔ_aff),
  convert(S.name.wrapper{T, 1}, dda.u_m_x_αΔ_aff),
  convert(S.name.wrapper{T, 1}, dda.s_l_αΔ_aff),
  convert(S.name.wrapper{T, 1}, dda.s_u_αΔ_aff),
  convert(S.name.wrapper{T, 1}, dda.rxs_l),
  convert(S.name.wrapper{T, 1}, dda.rxs_u),
)

function update_pt_aff!(
  x_m_l_αΔ_aff,
  u_m_x_αΔ_aff,
  s_l_αΔ_aff,
  s_u_αΔ_aff,
  Δxy_aff,
  Δs_l_aff,
  Δs_u_aff,
  x_m_lvar,
  uvar_m_x,
  s_l,
  s_u,
  α_aff_pri,
  α_aff_dual,
  ilow,
  iupp,
)
  x_m_l_αΔ_aff .= @views x_m_lvar .+ α_aff_pri .* Δxy_aff[ilow]
  u_m_x_αΔ_aff .= @views uvar_m_x .- α_aff_pri .* Δxy_aff[iupp]
  s_l_αΔ_aff .= s_l .+ α_aff_dual .* Δs_l_aff
  s_u_αΔ_aff .= s_u .+ α_aff_dual .* Δs_u_aff
end

# Mehrotra's Predictor-Corrector algorithm
function update_dd!(
  dda::DescentDirectionAllocsPC{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  pad::PreallocatedData{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}

  # solve system aff
  dda.Δxy_aff[1:(id.nvar)] .= .-res.rc
  dda.Δxy_aff[(id.nvar + 1):(id.nvar + id.ncon)] .= .-res.rb
  if typeof(pad) <: PreallocatedDataAugmented || typeof(pad) <: PreallocatedDataNormal
    dda.Δxy_aff[id.ilow] .+= pt.s_l
    dda.Δxy_aff[id.iupp] .-= pt.s_u
  elseif typeof(pad) <: PreallocatedDataNewton
    dda.Δs_l_aff .= .-itd.x_m_lvar .* pt.s_l
    dda.Δs_u_aff .= .-itd.uvar_m_x .* pt.s_u
  end

  cnts.w.write == true && write_system(cnts.w, pad.K, dda.Δxy_aff, :aff, cnts.k)
  out = solver!(dda.Δxy_aff, pad, dda, pt, itd, fd, id, res, cnts, T0, :aff)
  out == 1 && return out
  if typeof(pad) <: PreallocatedDataAugmented || typeof(pad) <: PreallocatedDataNormal
    dda.Δs_l_aff .= @views .-pt.s_l .- pt.s_l .* dda.Δxy_aff[id.ilow] ./ itd.x_m_lvar
    dda.Δs_u_aff .= @views .-pt.s_u .+ pt.s_u .* dda.Δxy_aff[id.iupp] ./ itd.uvar_m_x
  end

  if typeof(pt.x) <: Vector
    α_aff_pri, α_aff_dual = compute_αs(
      pt.x,
      pt.s_l,
      pt.s_u,
      fd.lvar,
      fd.uvar,
      dda.Δxy_aff,
      dda.Δs_l_aff,
      dda.Δs_u_aff,
      id.nvar,
    )
  else
    α_aff_pri, α_aff_dual = compute_αs_gpu(
      pt.x,
      pt.s_l,
      pt.s_u,
      fd.lvar,
      fd.uvar,
      dda.Δxy_aff,
      dda.Δs_l_aff,
      dda.Δs_u_aff,
      id.nvar,
      itd.store_vpri,
      itd.store_vdual_l,
      itd.store_vdual_u,
    )
  end

  # (x-lvar, uvar-x, s_l, s_u) .+= α_aff * Δ_aff
  update_pt_aff!(
    dda.x_m_l_αΔ_aff,
    dda.u_m_x_αΔ_aff,
    dda.s_l_αΔ_aff,
    dda.s_u_αΔ_aff,
    dda.Δxy_aff,
    dda.Δs_l_aff,
    dda.Δs_u_aff,
    itd.x_m_lvar,
    itd.uvar_m_x,
    pt.s_l,
    pt.s_u,
    α_aff_pri,
    α_aff_dual,
    id.ilow,
    id.iupp,
  )

  μ_aff =
    compute_μ(dda.x_m_l_αΔ_aff, dda.u_m_x_αΔ_aff, dda.s_l_αΔ_aff, dda.s_u_αΔ_aff, id.nlow, id.nupp)
  σ = (μ_aff / itd.μ)^3

  # corrector-centering step
  if typeof(pad) <: PreallocatedDataAugmented || typeof(pad) <: PreallocatedDataNormal
    dda.rxs_l .= @views -σ * itd.μ .+ dda.Δxy_aff[id.ilow] .* dda.Δs_l_aff
    dda.rxs_u .= @views σ * itd.μ .+ dda.Δxy_aff[id.iupp] .* dda.Δs_u_aff
    itd.Δxy .= 0
    itd.Δxy[id.ilow] .+= dda.rxs_l ./ itd.x_m_lvar
    itd.Δxy[id.iupp] .+= dda.rxs_u ./ itd.uvar_m_x
  elseif typeof(pad) <: PreallocatedDataNewton
    itd.Δxy[1:end] .= 0
    itd.Δs_l .= @views σ * itd.μ .- dda.Δxy_aff[id.ilow] .* dda.Δs_l_aff
    itd.Δs_u .= @views σ * itd.μ .+ dda.Δxy_aff[id.iupp] .* dda.Δs_u_aff
  end

  cnts.w.write == true && write_system(cnts.w, pad.K, itd.Δxy, :cc, cnts.k)
  out = solver!(itd.Δxy, pad, dda, pt, itd, fd, id, res, cnts, T0, :cc)
  out == 1 && return out
  if typeof(pad) <: PreallocatedDataAugmented || typeof(pad) <: PreallocatedDataNormal
    itd.Δs_l .= @views .-(dda.rxs_l .+ pt.s_l .* itd.Δxy[id.ilow]) ./ itd.x_m_lvar
    itd.Δs_u .= @views (dda.rxs_u .+ pt.s_u .* itd.Δxy[id.iupp]) ./ itd.uvar_m_x
  end

  # final direction
  itd.Δxy .+= dda.Δxy_aff
  itd.Δs_l .+= dda.Δs_l_aff
  itd.Δs_u .+= dda.Δs_u_aff

  return out
end

mutable struct DescentDirectionAllocsIPF{T <: Real, S} <: DescentDirectionAllocs{T, S}
  compl::S # complementarity s_lᵀ(x-lvar) + s_uᵀ(uvar-x)
  function DescentDirectionAllocsIPF(compl::AbstractVector{T}) where {T <: Real}
    S = typeof(compl)
    return new{T, S}(compl)
  end
end

DescentDirectionAllocsIPF(id::QM_IntData, S::DataType) =
  DescentDirectionAllocsIPF(S(undef, id.nvar))

convert(
  ::Type{<:DescentDirectionAllocs{T, S}},
  dda::DescentDirectionAllocsIPF{T0, S0},
) where {T <: Real, S, T0 <: Real, S0} =
  DescentDirectionAllocsIPF(convert(S.name.wrapper{T, 1}, dda.compl))

function update_dd!(
  dda::DescentDirectionAllocsIPF{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  pad::PreallocatedData{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}
  r, γ = T(0.999), T(0.05)
  # D = [s_l (x-lvar) + s_u (uvar-x)]
  dda.compl .= 0
  dda.compl[id.ilow] .+= pt.s_l .* itd.x_m_lvar
  dda.compl[id.iupp] .+= pt.s_u .* itd.uvar_m_x
  dda.compl[id.ifree] .= one(T)
  ξ = minimum(dda.compl) / itd.μ
  σ = γ * min((one(T) - r) * (one(T) - ξ) / ξ, T(2))^3

  itd.Δxy[1:(id.nvar)] .= .-res.rc
  itd.Δxy[(id.nvar + 1):(id.nvar + id.ncon)] .= .-res.rb
  if typeof(pad) <: PreallocatedDataAugmented || typeof(pad) <: PreallocatedDataNormal
    itd.Δxy[id.ilow] .+= pt.s_l - σ * itd.μ ./ itd.x_m_lvar
    itd.Δxy[id.iupp] .-= pt.s_u - σ * itd.μ ./ itd.uvar_m_x
  elseif typeof(pad) <: PreallocatedDataNewton
    itd.Δs_l .= σ * itd.μ .- itd.x_m_lvar .* pt.s_l
    itd.Δs_u .= σ * itd.μ .- itd.uvar_m_x .* pt.s_u
  end

  cnts.w.write == true && write_system(cnts.w, pad.K, itd.Δxy, :IPF, cnts.k)
  out = solver!(itd.Δxy, pad, dda, pt, itd, fd, id, res, cnts, T0, :IPF)
  out == 1 && return out
  if typeof(pad) <: PreallocatedDataAugmented || typeof(pad) <: PreallocatedDataNormal
    itd.Δs_l .= @views (σ * itd.μ .- pt.s_l .* itd.Δxy[id.ilow]) ./ itd.x_m_lvar .- pt.s_l
    itd.Δs_u .= @views (σ * itd.μ .+ pt.s_u .* itd.Δxy[id.iupp]) ./ itd.uvar_m_x .- pt.s_u
  end
end


text:
alphaOnlyMatchRegex = r"^[a-zA-Z0-9_]*$"

"""
$(SIGNATURES)
Returns the transaction for a given query.
NOTE: Must commit(transaction) after you're done.
"""
function _queryNeo4j(neo4jInstance::Neo4jInstance, query::String)
    loadtx = transaction(neo4jInstance.connection)
    result = loadtx(query; submit=true)
    if length(result.errors) > 0
        error(string(result.errors))
    end
    # Have to finish the transaction
    commit(loadtx)
    return result
end

"""
$(SIGNATURES)
Create a node and optionally specify a parent.
Note: Using symbols so that the labels obey Neo4j requirements
"""
function _createNode(neo4jInstance::Neo4jInstance, labels::Vector{String}, properties::Dict{String, Any}, parentNode::Union{Nothing, Neo4j.Node}, relationshipLabel::Symbol=:NOTHING)::Neo4j.Node
    createdNode = Neo4j.createnode(neo4jInstance.graph, properties)
    addnodelabels(createdNode, labels)
    parentNode == nothing && return createdNode
    # Otherwise create the relationship
    createrel(parentNode, createdNode, String(relationshipLabel))
    return createdNode
end

"""
$(SIGNATURES)
Returns the list of CloudGraph nodes that matches the Cyphon query.
The nodes of interest should be labelled 'node' because the query will use the return of id(node)
#Example
matchCondition = '(u:User)-[:ROBOT]->(n:Robot:NewRobot)-[:SESSION]->(node:Session)'.
Will return all sessions because of 'node:Session'.
If orderProperty is not ==, then 'order by n.{orderProperty} will be appended to the query'.
So can make orderProperty = label or id.
"""
function _getLabelsFromCyphonQuery(neo4jInstance::Neo4jInstance, matchCondition::String, orderProperty::String="")::Vector{Symbol}
    query = "match $matchCondition return distinct(node.label) $(orderProperty != "" ? "order by node.$orderProperty" : "")";
    @debug "[Query] $query"
    result = _queryNeo4j(neo4jInstance, query)
    nodeIds = map(node -> node["row"][1], result.results[1]["data"])
    return Symbol.(nodeIds)
end

"""
$(SIGNATURES)
Get a node property - returns nothing if not found
"""
function _getNodeProperty(neo4jInstance::Neo4jInstance, nodeLabels::Vector{String}, property::String)
    query = "match (n:$(join(nodeLabels, ":"))) return n.$property"
    @debug "[Query] $query"
    result = _queryNeo4j(neo4jInstance, query)
    length(result.results[1]["data"]) != 1 && error("No data returned from the query.")
    length(result.results[1]["data"][1]["row"]) != 1 && error("No data returned from the query.")
    return result.results[1]["data"][1]["row"][1]
end

"""
$(SIGNATURES)
Set a node property - returns count of changed nodes.
"""
function _setNodeProperty(neo4jInstance::Neo4jInstance, nodeLabels::Vector{String}, property::String, value::Union{String, Int, Float64})
    if value isa String
        value = "\""*replace(value, "\"" => "\\\"")*"\"" # Escape strings
    end
    query = """
    match (n:$(join(nodeLabels, ":")))
    set n.$property = $value
    return count(n)"""
    @debug "[Query] $query"
    result = _queryNeo4j(neo4jInstance, query)
    length(result.results[1]["data"]) != 1 && return 0
    length(result.results[1]["data"][1]["row"]) != 1 && return 0
    return result.results[1]["data"][1]["row"][1]
end

"""
$(SIGNATURES)
Get a node's tags
"""
function _getNodeTags(neo4jInstance::Neo4jInstance, nodeLabels::Vector{String})::Union{Nothing, Vector{String}}
    query = "match (n:$(join(nodeLabels, ":"))) return labels(n)"
    result = _queryNeo4j(neo4jInstance, query)
    length(result.results[1]["data"]) != 1 && return nothing
    return result.results[1]["data"][1]["row"][1]
end

function _getNodeCount(neo4jInstance::Neo4jInstance, nodeLabels::Vector{String})::Int
    query = "match (n:$(join(nodeLabels, ":"))) return count(n)"
    result = _queryNeo4j(neo4jInstance, query)
    length(result.results[1]["data"]) != 1 && return 0
    length(result.results[1]["data"][1]["row"]) != 1 && return 0
    return result.results[1]["data"][1]["row"][1]
end

"""
$(SIGNATURES)
Returns the list of CloudGraph nodes that matches the Cyphon query.
The nodes of interest should be labelled 'node' because the query will use the return of id(node)
#Example
matchCondition = '(u:User)-[:ROBOT]->(n:Robot:NewRobot)-[:SESSION]->(node:Session)'.
Will return all sessions because of 'node:Session'.
If orderProperty is not ==, then 'order by n.{orderProperty} will be appended to the query'.
So can make orderProperty = label or id.
"""
function _getNeoNodesFromCyphonQuery(neo4jInstance::Neo4jInstance, matchCondition::String, orderProperty::String="")::Vector{Neo4j.Node}
    # 2. Perform the transaction
    loadtx = transaction(neo4jInstance.connection)
    query = "match $matchCondition return id(node) $(orderProperty != "" ? "order by node.$orderProperty" : "")";
    nodes = loadtx(query; submit=true)
    if length(nodes.errors) > 0
        error(string(nodes.errors))
    end
    # 3. Format the result
    nodes = map(node -> getnode(neo4jInstance.graph, node["row"][1]), nodes.results[1]["data"])
    # Have to finish the transaction
    commit(loadtx)
    return nodes
end

"""
$(SIGNATURES)
Utility function to get a Neo4j node sessions for a robot.
"""
function _getSessionNeoNodesForRobot(neo4jInstance::Neo4jInstance, userId::String, robotId::String)::Vector{Neo4j.Node}
    # 2. Perform the transaction
    nodes = _getNeoNodesFromCyphonQuery(neo4jInstance, "(u:USER:$userId)-[:ROBOT]->(r:ROBOT:$robotId)-[:SESSION]->(node:SESSION)", "id")
    return nodes
end

"""
$(SIGNATURES)
Bind the SESSION node to the inital variable.
Checks for existence.
"""
function _bindSessionNodeToInitialVariable(neo4jInstance::Neo4jInstance, userId::String, robotId::String, sessionId::String, initialVariableLabel::String)::Nothing
    # 2. Perform the transaction
    loadtx = transaction(neo4jInstance.connection)
    query = """
    match (session:SESSION:$userId:$robotId:$sessionId),
    (var:VARIABLE:$userId:$robotId:$sessionId
    {label: '$initialVariableLabel'})
    WHERE NOT (session)-[:VARIABLE]->()
    CREATE (session)-[:VARIABLE]->(var) return id(var)
    """;
    loadtx(query; submit=true)
    commit(loadtx)
    return nothing
end

"""
$(SIGNATURES)
Try get a Neo4j node ID from a node label.
"""
function _tryGetNeoNodeIdFromNodeLabel(neo4jInstance::Neo4jInstance, userId::String, robotId::String, sessionId::String, nodeLabel::Symbol; nodeType::Union{Nothing, String}=nothing)::Union{Nothing, Int}
    @debug "Looking up symbolic node ID where n.label = '$nodeLabel'..."
    if nodeType == nothing
        nodes = _getNeoNodesFromCyphonQuery(neo4jInstance, "(node:$userId:$robotId:$sessionId) where exists(node.label) and node.label = \"$(string(nodeLabel))\"")
    else
        nodes = _getNeoNodesFromCyphonQuery(neo4jInstance, "(node:$nodeType:$userId:$robotId:$sessionId) where exists(node.label) and node.label = \"$(string(nodeLabel))\"")
    end

    if length(nodes) != 1
        return nothing
    end
    nodeId = nodes[1].id
    @debug "Found NeoNode ID = $nodeId..."
    return nodeId
end

"""
$(SIGNATURES)
Gets the searchable labels for any CGDFG type.
"""

"""
$(SIGNATURES)

Build a Cypher-compliant set of properies from a JSON dictionary.
Note individual values are serialized if they are not already.
"""
function _dictToNeo4jProps(dict::Dict{String, Any})::String
    # TODO: Use an IO buffer/stream for this
    return "{" * join(map((k) -> "$k: $(JSON.json(dict[k]))", collect(keys(dict))), ", ")*"}"
end

"""
$(SIGNATURES)

Get the Neo4j labels for any node type.
"""
function _getLabelsForType(dfg::CloudGraphsDFG, type::Type; parentKey::Union{Nothing, Symbol}=nothing)::Vector{String}
    # Simple validation
    isempty(dfg.userId) && error("The DFG object's userID is empty, please specify a user ID.")
    isempty(dfg.robotId) && error("The DFG object's robotID is empty, please specify a robot ID.")
    isempty(dfg.sessionId) && error("The DFG object's sessionId is empty, please specify a session ID.")

    labels = []
    type == User && (labels = [dfg.userId, "USER"])
    type == Robot && (labels = [dfg.userId, dfg.robotId, "ROBOT"])
    type == Session && (labels = [dfg.userId, dfg.robotId, dfg.sessionId, "SESSION"])
    type <: DFGVariable &&
        (labels = [dfg.userId, dfg.robotId, dfg.sessionId, "VARIABLE"])
    type <: DFGFactor &&
        (labels = [dfg.userId, dfg.robotId, dfg.sessionId, "FACTOR"])
    type <: AbstractPointParametricEst &&
        (labels = [dfg.userId, dfg.robotId, dfg.sessionId, "PPE"])
    type <: AbstractBigDataEntry &&
        (labels = [dfg.userId, dfg.robotId, dfg.sessionId, "DATA"])
    # Some are children of nodes, so add that in if it's set.
    parentKey != nothing && push!(labels, String(parentKey))
    return labels
end

"""
$(SIGNATURES)

Get the Neo4j labels for any node instance.
"""
function _getLabelsForInst(dfg::CloudGraphsDFG,
                            inst::Union{User, Robot, Session, N, APPE, ABDE}; parentKey::Union{Nothing, Symbol}=nothing)::Vector{String} where
                            {N <: DFGNode, APPE <: AbstractPointParametricEst, ABDE <: AbstractBigDataEntry}
    labels = _getLabelsForType(dfg, typeof(inst), parentKey=parentKey)
    typeof(inst) <: DFGVariable && push!(labels, String(getLabel(inst)))
    typeof(inst) <: DFGFactor && push!(labels, String(getLabel(inst)))
    typeof(inst) <: AbstractPointParametricEst && push!(labels, String(inst.solverKey))
    typeof(inst) <: AbstractBigDataEntry && push!(labels, String(inst.key))
    # Some are children of nodes, so add that in if it's set.
    return labels
end


text:
function EcksteinHechlerPropagatorBuilder(arg0::Orbit, arg1::UnnormalizedSphericalHarmonicsProvider, arg2::PositionAngle, arg3::jdouble)
    return EcksteinHechlerPropagatorBuilder((Orbit, UnnormalizedSphericalHarmonicsProvider, PositionAngle, jdouble), arg0, arg1, arg2, arg3)
end

function EcksteinHechlerPropagatorBuilder(arg0::Orbit, arg1::UnnormalizedSphericalHarmonicsProvider, arg2::PositionAngle, arg3::jdouble, arg4::AttitudeProvider)
    return EcksteinHechlerPropagatorBuilder((Orbit, UnnormalizedSphericalHarmonicsProvider, PositionAngle, jdouble, AttitudeProvider), arg0, arg1, arg2, arg3, arg4)
end

function EcksteinHechlerPropagatorBuilder(arg0::Orbit, arg1::jdouble, arg10::PositionAngle, arg11::jdouble, arg2::jdouble, arg3::TideSystem, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::OrbitType)
    return EcksteinHechlerPropagatorBuilder((Orbit, jdouble, PositionAngle, jdouble, jdouble, TideSystem, jdouble, jdouble, jdouble, jdouble, jdouble, OrbitType), arg0, arg1, arg10, arg11, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function build_propagator(obj::EcksteinHechlerPropagatorBuilder, arg0::Vector{jdouble})
    return jcall(obj, "buildPropagator", Propagator, (Vector{jdouble},), arg0)
end


text:
__precompile__()

module PositiveFactorizations

using Compat
using Compat.view

export Positive

immutable Positive{T<:Real} end

include("cholesky.jl")
include("eig.jl")

end # module


text:
module Oolong

const OL = Oolong
export OL

include("core.jl")
include("parameter_server.jl")
include("serve.jl")

function __init__()
    init()
end

end


text:
module Digitization
using ParallelKMeans
using Statistics

export digitization_k, Model

mutable struct Model
   centers::Array{Float64,}
   hashmap::Dict{Char,Int64}
end

function digitization_k(array, k::Int, nthreads::Int)
    features = collect(Matrix(array[:, 1:2])');
    result = kmeans(features, k; n_threads=nthreads);
    labels = result.assignments;
    symbols = map(x -> Char(x+33), labels);
    hashmap = Dict(zip(symbols, labels));
    centers = zeros((0, 2));
    features = features';
    for i in sort(unique(labels))
        condition = (labels .== i);
        group = features[condition, :];
        center = mean!([1. 1.], group);
        centers = vcat(centers, center);
    end
    model = Model(centers, hashmap);
    return symbols, model;
end

end

max_stars_repo_name: nla-group/Julia-ABBA

max_issues_repo_name: nla-group/Julia-ABBA

max_forks_repo_name: nla-group/Julia-ABBA


text:
function compute_JY(z, n, dim)
    if dim == 2
        J = besselj.(n, z)
        Y = bessely.(n, z)
        dJ = (besselj.(n - 1, z) - besselj.(n + 1, z)) / 2
        dY = (bessely.(n - 1, z) - bessely.(n + 1, z)) / 2
    else
        J = sphericalbesselj.(n, z)
        Y = sphericalbessely.(n, z)

        dJ = [copy(z);]
        dY = [copy(z);]
        ind = z .> 0
        zz = z[ind]
        ind = [ind;]
        dJ[ind] .=
            (
                sphericalbesselj.(n - 1, zz) - sphericalbesselj.(n + 1, zz) -
                sphericalbesselj.(n, zz) ./ zz
            ) / 2
        dY[ind] .=
            (
                sphericalbessely.(n - 1, zz) - sphericalbessely.(n + 1, zz) -
                sphericalbessely.(n, zz) ./ zz
            ) / 2

        ind = [z .== 0;]
        dJ[ind] .= (n == 1) / 3
        dY[ind] .= Inf

        if size(z) == ()
            J, Y, dJ, dY = J[1], Y[1], dJ[1], dY[1]
        end
    end

    J, Y, dJ, dY
end


text:
# coding:utf-8

doc=
"""
  レファレンス質問収集プログラム
  レファレンス協同データベースより
"""

using DataFrames
using Requests

# 検索結果からidをパース
function id_parse(qid, html, result_base)
  for m in matchall(result_base, html)
    id = match(result_base, m)[1]
    println(id)
    push!(qid, id)
  end
  return qid
end

# collect_id
function collect_id(id_list)
  # base url
  search_base ="http://crd.ndl.go.jp/reference/modules/d3ndlcrdsearch/index.php?page=detail_list&type=reference&mcmd=200&st=update&asc=desc&dtltbs=1&ndc_lk=1&ndc1="
  # search result
  result_base = r"http://crd.ndl.go.jp/reference/modules/d3ndlcrdentry/index.php\?page=ref_view\&ldtl=1\&dtltbs=\d\&mcmd=200\&st=update\&asc=desc(?:|\&pg=\d+)\&ndc1=\d\&ndc_lk=1\&id=(\d+)"
  # NDC一次区分
  ndc_top = collect(0:9)
  qid = []
  # question id
  for ndc in ndc_top
    println(ndc)
    pg = 1
    while true
      println(search_base*"$ndc"*"&pg=$pg")
      html = readall(get(search_base*"$ndc"*"&pg=$pg"))
      qid, pg = id_parse(qid, html, result_base), pg + 1
      !contains(html, "&pg=$pg") && break
    end
  end
  gc()
  f = open(id_list, "w")
  write(f, join(qid, "\n"))
  close(f)
end

# レファレンス記録の取得&パース
function get_reference_question(id_list, page_dict)
  # 保存用
  !isdir("jp_reference") && mkdir("jp_reference")
  # 確認用
  if isfile(page_dict)
    d = open(deserialize, page_dict)
  else
    d = Dict()
  end
  # base URL
  base = "http://crd.ndl.go.jp/reference/detail?page=ref_view&id="
  # レファレンス記録取得
  qid = unique(split(readall(id_list), "\n"))
  try
    for (i, id) in enumerate(qid)
      println(i,": ",id)
      haskey(d, id) && continue
      html = readall(get(base*"$id"))
      m = match(r"
(.+?)
", html)[1]
      question = lstrip(replace(m, r"<.+?>|\n", ""))
      d[id] = [question, "jp_reference/$id.html"]
      println(question)
      write(open("jp_reference/$id.html", "w"), html)
      i % 100 == 0 && open(io -> serialize(io, d), page_dict, "w")
      sleep(rand()*0.9)
    end
  catch
    # エラーハンドリン
    println("error?")
　  open(io -> serialize(io, d), page_dict, "w")
    # 再送
    println("sleep ...")
    sleep(10)
    get_reference_question(id_list, page_dict)
  end
end

# page_idを取得
id_list = "jp_reference_id.dat"
collect_id(id_list)

# 記事取得
page_dict = "jp_reference_description.dict"
get_reference_question(id_list, page_dict)

max_stars_repo_name: kwnsiy/ReferenceQuestionQrawaler

max_issues_repo_name: kwnsiy/ReferenceQuestionQrawaler

max_forks_repo_name: kwnsiy/ReferenceQuestionQrawaler


text:
#include("HPDE_options.jl")
##########################################################

mutable struct FVMPDEProblem
    grid::FVMPDEGrid

    """
    Flux function F in x direction

    """
    F::Function

    """
    Flux function G in y direction
    !!! compat
        Will be added in future versions
    """
    G::Function

    """
    Flux function H in z direction
    !!! compat
        Will be added in future versions
    """
    H::Function

    """
    Source function
    !!! compat
        Will be added in future versions
    """
    S::Function

    nvars::Int8

    param::Any

    U0

    U_min
    U_max

    jac
    ode_problem_kwargs
    tspan
    t

    function FVMPDEProblem(grid::FVMPDEGrid, tspan, U0; kwargs...)
        if haskey(kwargs,:jacobian)
            jac = kwargs[:jacobian]
            kwargs = Dict([p for p in pairs(kwargs) if p[1] != :jacobian])
        else
            jac = true
        end
        zeroFunc(U) = zeros(size(U))
        F = kwargs[:F]
        @assert size(U0,1) == grid.nx
        if haskey(kwargs,:G)
            G = kwargs[:G]
            kwargs = Dict([p for p in pairs(kwargs) if p[1] != :G])
            @assert size(U0,2) == grid.ny
            if haskey(kwargs,:H)
                H = kwargs[:H]
                kwargs = Dict([p for p in pairs(kwargs) if p[1] != :H])
                nvars = size(U0,4)
                @assert grid.dimension == _3D
                @assert length(size(U0))-1 == 3
                @assert size(U0,3) == grid.nz
            else
                H = zeroFunc
                nvars = size(U0,3)
                @assert grid.dimension == _2D
                @assert length(size(U0))-1 == 2
            end
        else
            @assert !haskey(kwargs,:H)
            G = zeroFunc
            H = zeroFunc
            nvars = size(U0,2)
            @assert grid.dimension == _1D
            @assert length(size(U0))-1 == 1
        end

        if haskey(kwargs,:source)
            S = kwargs[:source]
            kwargs = Dict([p for p in pairs(kwargs) if p[1] != :source])
        else
            S = zeroFunc
        end
        if jac
            @assert size(F(U0[grid.indices[1],:])) == size(U0[grid.indices[1],:])
            @assert size(G(U0[grid.indices[1],:])) == size(U0[grid.indices[1],:])
            @assert size(H(U0[grid.indices[1],:])) == size(U0[grid.indices[1],:])
        else
            @assert size(F(U0[grid.indices[1],:]),1) == length(U0[grid.indices[1],:])
            @assert size(F(U0[grid.indices[1],:]),2) == length(U0[grid.indices[1],:])
            # TODO: add @assert for G and H
        end
        @assert size(S(U0)) == size(U0)

        #grid.U = reshape(U0,grid.nx,grid.ny,grid.nz,nvars)
        #grid.U = copy(U0)
        if haskey(kwargs,:parameters)
            θ = kwargs[:parameters]
            kwargs = Dict([p for p in pairs(kwargs) if p[1] != :parameters])
        else
            θ = nothing
        end
        if haskey(kwargs,:ul)
            U_min = kwargs[:ul]
            kwargs = Dict([p for p in pairs(kwargs) if p[1] != :ul])
        else
            U_min = -Inf * ones(size(U0))
        end
        if haskey(kwargs,:uh)
            U_max = kwargs[:uh]
            kwargs = Dict([p for p in pairs(kwargs) if p[1] != :uh])
        else
            U_max = Inf * ones(size(U0))
        end

        @assert size(U_min) == size(U0)
        @assert size(U_max) == size(U0)

        #@assert checkDim(D) == size(U0,)

        return new(grid,F, G, H, S, nvars, θ, U0, U_min, U_max,jac,kwargs,tspan,0)
    end

end


text:
using BenchmarkTools, NDTools

function compare_reverse(arr, dims=Tuple(1:ndims(arr)))
    print("## runtime with collect\nsize(arr) = $(size(arr))\n")
    print("\treverse:\t")
    @btime reverse($arr, dims=$dims)
    print("\treverse_view:\t")
    @btime collect(reverse_view($arr, dims=$dims))
    print("\n")
end

function compare_reverse_broadcast(arr, dims=Tuple(1:ndims(arr)))
    print("## Broadcast \nsize(arr) = $(size(arr))\n")
    print("\treverse:\t")
    arr2 = copy(arr)
    @btime $arr2 .= abs2.(reverse!($arr, dims=$dims))
    print("\treverse_view:\t")
    @btime $arr2 .= abs2.(reverse_view($arr, dims=$dims))
    print("\n")
end

function main()
    compare_reverse_broadcast(randn(10, 10))
    compare_reverse_broadcast(randn(1013, 1024))
    compare_reverse_broadcast(randn(512, 13))
    compare_reverse_broadcast(randn(512, 133, 231,2))
    compare_reverse_broadcast(randn(512, 133, 231,2), (1,4))

    compare_reverse(randn(10, 10))
    compare_reverse(randn(1013, 1024))
    compare_reverse(randn(512, 13))
    compare_reverse(randn(512, 133, 231,2))
    compare_reverse(randn(512, 133, 231,2), (1,4))
end


text:
using Test, WAVI

@testset "WAVI tests" begin
    @testset "Pos Fraction" begin
        @info "testing pos fraction"
        z=[-1.0 -1.0 -1.0;-1.0 1.0 -1.0;-1.0 -1.0 -1.0]
        pfh,pfu,pfv=WAVI.pos_fraction(z)
        @test pfh == [0.0 0.0 0.0;0.0 0.5 0.0;0.0 0.0 0.0]
        @test pfu == [0.0 0.0 0.0;0.0 0.25 0.0;0.0 0.25 0.0;0.0 0.0 0.0]
        @test pfv == [0.0 0.0 0.0 0.0;0.0 0.25 0.25 0.0;0.0 0.0 0.0 0.0]
    end

    @testset "Pos Fraction Mask" begin
        @info "testing post fraction mask"
        z=[-1.0 -1.0 -1.0;-1.0 1.0 -1.0;-1.0 -1.0 -1.0]
        mask=[false false false; false true false; false false false]
        pfh,pfu,pfv=WAVI.pos_fraction(z,mask=mask)
        @test pfh == [0.0 0.0 0.0;0.0 1.0 0.0;0.0 0.0 0.0]
        @test pfu == [0.0 0.0 0.0;0.0 0.5 0.0;0.0 0.5 0.0;0.0 0.0 0.0]
        @test pfv == [0.0 0.0 0.0 0.0;0.0 0.5 0.5 0.0;0.0 0.0 0.0 0.0]
    end

    @testset "Height and volume above floatation" begin
        @info "testing height and volume above floatation functions"
        grid = Grid(nx = 20, ny = 25, dx = 1000., dy = 1000.)
        bed = -200.0 * ones(grid.nx, grid.ny)
        initial_conditions = InitialConditions(initial_thickness = 500.0*ones(grid.nx, grid.ny))
        model = Model(grid = grid, bed_elevation = bed, initial_conditions = initial_conditions)
        update_state!(model)
        @test volume_above_floatation(model.fields.gh.h, model.fields.gh.b, Ref(model.params), model.grid) ≈ (500.0 - (1028.0/918.0)*(200.0) ) .* 25 .* 20 .* 1000.0 .^2

    end
end