text:
#importall CUTEst
#using SparseArrays
#import CUTEst
include("../parameters.jl")
include("../linear_system_solvers/linear_system_solvers.jl")
include("../kkt_system_solver/include.jl")
include("../line_search/line_search.jl")
include("../IPM/ipm.jl")
include("../JuMPinterface.jl")
export get_original_x, get_y, number_constraints

mutable struct Class_bounds
    l_i::Array{Int64,1}
    u_i::Array{Int64,1}
    l::Array{Float64,1}
    u::Array{Float64,1}

    function Class_bounds(lb_vec::Array{Float64,1}, ub_vec::Array{Float64,1})
        @assert(length(lb_vec) == length(ub_vec))
        this = new([],[],[],[]);
        for i = 1:length(lb_vec)
            if lb_vec[i] > -Inf
               push!(this.l_i, i)
               push!(this.l, lb_vec[i])
            end

            if ub_vec[i] < Inf
               push!(this.u_i, i)
               push!(this.u, ub_vec[i])
            end
        end

        return this;
    end
end

function _i_not_fixed(m::NLPModels.AbstractNLPModel)
    return (1:m.meta.nvar)[m.meta.lvar .!= m.meta.uvar]
end

function _i_not_fixed(variable_info::Vector{VariableInfo})
    vector_indixes = zeros(Int64, 0)
    for i in 1:length(variable_info)
        if !variable_info[i].is_fixed
            push!(vector_indixes, i)
        end
    end
    return vector_indixes
end

function extract_lcon_ucon(solver::OnePhaseSolver)
    lcon = zeros(Float64, 0)
    ucon = zeros(Float64, 0)

    #Linear constraints
    for i in 1:length(solver.linear_eq_constraints)
        lower = solver.linear_eq_constraints[i].set.value
        push!(lcon, lower)
        push!(ucon, lower)
    end

    for i in 1:length(solver.linear_ge_constraints)
        lower = solver.linear_ge_constraints[i].set.lower
        upper = Inf
        push!(lcon, lower)
        push!(ucon, upper)
    end

    for i in 1:length(solver.linear_le_constraints)
        lower = -Inf
        upper = solver.linear_le_constraints[i].set.upper
        push!(lcon, lower)
        push!(ucon, upper)
    end

    for i in 1:length(solver.linear_int_constraints)
        lower = solver.linear_int_constraints[i].set.lower
        upper = solver.linear_int_constraints[i].set.upper
        push!(lcon, lower)
        push!(ucon, upper)
    end

    #Quadratic Constraints
    for i in 1:length(solver.quadratic_eq_constraints)
        lower = solver.quadratic_eq_constraints[i].set.lower
        push!(lcon, lower)
        push!(ucon, lower)
    end

    for i in 1:length(solver.quadratic_ge_constraints)
        lower = solver.quadratic_ge_constraints[i].set.lower
        upper = Inf
        push!(lcon, lower)
        push!(ucon, upper)
    end

    for i in 1:length(solver.quadratic_le_constraints)
        lower = -Inf
        upper = solver.quadratic_le_constraints[i].set.upper
        push!(lcon, lower)
        push!(ucon, upper)
    end

    for i in 1:length(solver.quadratic_int_constraints)
        lower = solver.quadratic_int_constraints[i].set.lower
        upper = solver.quadratic_int_constraints[i].set.upper
        push!(lcon, lower)
        push!(ucon, upper)
    end

    #Nonlinear Constraints
    for i in 1:length(solver.nlp_data.constraint_bounds)
        push!(lcon, solver.nlp_data.constraint_bounds[i].lower)
        push!(ucon, solver.nlp_data.constraint_bounds[i].upper)
    end

    return lcon, ucon
end

function extract_lvar_uvar(variable_info::Vector{VariableInfo}, ind::Vector{Int64})
    lvar = zeros(Float64, 0)
    uvar = zeros(Float64, 0)
    for i in ind
        push!(lvar, variable_info[i].lower_bound)
        push!(uvar, variable_info[i].upper_bound)
    end
    return lvar, uvar
end

mutable struct Class_CUTEst <: abstract_nlp
    nlp::Union{NLPModels.AbstractNLPModel, Nothing}
    solver::Union{OnePhaseSolver, Nothing}

    bcon::Class_bounds
    bvar::Class_bounds

    function Class_CUTEst(nlp::NLPModels.AbstractNLPModel)
        ind = _i_not_fixed(nlp)
        return new(nlp, nothing, Class_bounds(nlp.meta.lcon[:], nlp.meta.ucon[:]), Class_bounds(nlp.meta.lvar[ind], nlp.meta.uvar[ind]))
    end

    function Class_CUTEst(solver::OnePhaseSolver)
        ind = _i_not_fixed(solver.variable_info)
        lcon, ucon = extract_lcon_ucon(solver)
        lvar, uvar = extract_lvar_uvar(solver.variable_info, ind)
        return new(nothing, solver, Class_bounds(lcon, ucon), Class_bounds(lvar, uvar))
    end
end

function suggested_starting_point(m::Class_CUTEst)
    if m.nlp != nothing
        ind = _i_not_fixed(m.nlp)
        return deepcopy(m.nlp.meta.x0[ind])
    end

    ind = _i_not_fixed(m.solver.variable_info)
    non_fixed_variables = m.solver.variable_info[ind]
    starting_points_vector = zeros(0)
    for variable in non_fixed_variables
        push!(starting_points_vector, variable.start == nothing ? 0.0 : variable.start)
    end

    return deepcopy(starting_points_vector)
end

function ncon(m::Class_CUTEst)
    return nbounds_orginal(m) + ncons_orginal(m)
end

function lb(x::Array{Float64,1}, bd::Class_bounds)
    return x[bd.l_i] - bd.l
end

function ub(x::Array{Float64,1}, bd::Class_bounds)
    return bd.u - x[bd.u_i]
end

function number_constraints(solver::OnePhaseSolver)
    linear_constraints_count = number_constraints_linear(solver)
    quadratic_constraints_count = number_constraints_quadratic(solver)
    non_linear_constraints_count = number_constraints_non_linear(solver)
    return linear_constraints_count + quadratic_constraints_count + non_linear_constraints_count
end

function number_constraints_linear(solver::OnePhaseSolver)
    less_than_constraints_count = length(solver.linear_le_constraints)
    greater_than_constraints_count = length(solver.linear_ge_constraints)
    equality_constraints_count = length(solver.linear_eq_constraints)
    interval_constraints_count = length(solver.linear_int_constraints)
    return less_than_constraints_count + greater_than_constraints_count + equality_constraints_count + interval_constraints_count
end

function number_constraints_quadratic(solver::OnePhaseSolver)
    less_than_constraints_count = length(solver.quadratic_le_constraints)
    greater_than_constraints_count = length(solver.quadratic_ge_constraints)
    equality_constraints_count = length(solver.quadratic_eq_constraints)
    interval_constraints_count = length(solver.quadratic_int_constraints)
    return less_than_constraints_count + greater_than_constraints_count + equality_constraints_count + interval_constraints_count
end

function number_constraints_non_linear(solver::OnePhaseSolver)
    return length(solver.nlp_data.constraint_bounds)
end

function linear_cons(m::Class_CUTEst)
    if m.nlp != nothing
        is_linear = zeros(m.nlp.meta.ncon)
        is_linear[m.nlp.meta.lin] .= 1.0
        vec = [is_linear[m.bcon.l_i]; is_linear[m.bcon.u_i]; ones(nbounds_orginal(m))]
        return 1.0 .== vec
    end

    is_linear = zeros(number_constraints(m.solver))
    is_linear[[i for i in 1:number_constraints_linear(m.solver)]] .= 1.0
    vec = [is_linear[m.bcon.l_i]; is_linear[m.bcon.u_i]; ones(nbounds_orginal(m))]
    return 1.0 .== vec
end

function ineq_cons(m::Class_CUTEst)
    if m.nlp != nothing
        is_ineq = zeros(m.nlp.meta.ncon)
        is_ineq[m.nlp.meta.lcon .== m.nlp.meta.ucon] .= 1.0
        vec = [is_ineq[m.bcon.l_i]; is_ineq[m.bcon.u_i]; ones(nbounds_orginal(m))]
        return 1.0 .== vec
    end

    solver = m.solver
    is_ineq =  zeros(number_constraints(solver))

    index = 1 + length(solver.linear_le_constraints) + length(solver.linear_ge_constraints)

    is_linear[[i for i in index:length(solver.linear_eq_constraints)]] .= 1.0

    index += length(solver.linear_eq_constraints)

    for i in 1:length(solver.linear_int_constraints)
        lower = solver.linear_int_constraints[i].set.lower
        upper = solver.linear_int_constraints[i].set.upper
        if lower == upper
            is_ineq[index] .= 1.0
        end
        index +=1
    end

    index += length(solver.quadratic_le_constraints) + length(solver.quadratic_ge_constraints)

    is_linear[[i for i in index:length(solver.quadratic_eq_constraints)]] .= 1.0

    index += length(solver.quadratic_eq_constraints)

    for i in 1:length(solver.quadratic_int_constraints)
        lower = solver.quadratic_int_constraints[i].set.lower
        upper = solver.quadratic_int_constraints[i].set.upper
        if lower .== upper
            is_ineq[index] .= 1.0
        end
        index +=1
    end

    #Nonlinear Constraints
    for i in 1:length(solver.nlp_data.constraint_bounds)
        lower = solver.nlp_data.constraint_bounds[i].lower
        upper = solver.nlp_data.constraint_bounds[i].upper
        if lower .== upper
            is_ineq[index] .= 1.0
        end
        index +=1
    end

    vec = [is_ineq[m.bcon.l_i]; is_ineq[m.bcon.u_i]; ones(nbounds_orginal(m))]
    return 1.0 .== vec
end

function nbounds_orginal(nlp::Class_CUTEst)
    return length(nlp.bvar.l_i) + length(nlp.bvar.u_i)
end

function ncons_orginal(nlp::Class_CUTEst)
    return length(nlp.bcon.l_i) + length(nlp.bcon.u_i)
end

function cons_indicies(nlp::Class_CUTEst)
    m = ncons_orginal(nlp)
    if m > 0
      return 1:m
    else
      return []
    end
end

function bound_indicies(nlp::Class_CUTEst)
    m = ncons_orginal(nlp)
    r = nbounds_orginal(nlp)
    if r > 0
      return (m + 1):(m + r)
    else
      return [];
    end
end

function y_l_con(y::Array{Float64,1}, m::Class_CUTEst)
    return y[1:length(m.bcon.l)]
end

function y_u_con(y::Array{Float64,1}, m::Class_CUTEst)
    n_lcon = length(m.bcon.l)
    return y[(n_lcon + 1):(n_lcon + length(m.bcon.u))]
end

function eval_f(m::Class_CUTEst, x::Array{Float64,1})
    if m.nlp != nothing
        return obj(m.nlp, _cute_x(m, x) )
    end

    return MOI.eval_objective(m.solver.nlp_data.evaluator, _cute_x(m, x))
end

function evaluateScalarAffineFunctionConstraints(constraints::Vector{OnePhase.ConstraintInfo{MathOptInterface.ScalarAffineFunction{Float64}}}, x::Array{Float64,1})
    constraint_values = zeros(Float64, 0)
    for constraint in constraints
        value = constraint.func.constant
        for term in constraint.func.terms
            value += term.coefficient * x[term.variable.value]
        end
        push!(constraint_values, value)
    end
    return constraint_values
end

function evaluateScalarQuadraticFunctionConstraints(constraints::Vector{OnePhase.ConstraintInfo{MathOptInterface.ScalarQuadraticFunction{Float64}}}, x::Array{Float64,1})
    constraint_values = zeros(Float64, 0)
    for constraint in constraints
        value = constraint.func.constant
        for term in constraint.func.affine_terms
            value += term.coefficient * x[term.variable.value]
        end
        for term in constraint.func.quadratic_terms
            value += term.coefficient * x[term.variable_1.value] * x[term.variable_2.value]
        end
        push!(constraint_values, value)
    end
    return constraint_values
end

function evaluateScalarAffineFunctionConstraintsJacobianStructure(constraints::Vector{OnePhase.ConstraintInfo{MathOptInterface.ScalarAffineFunction{Float64}}}, x::Array{Float64,1})
    jac_structures = nothing
    count = 1
    for constraint in constraints
        jac_structure = zeros(Float64, length(x))
        for term in constraint.func.terms
            index = term.variable.value
            value = term.coefficient
            jac_structure[index] = value
        end
        if count == 1
            jac_structures = jac_structure
        else
            jac_structures = hcat(jac_structures, jac_structure)
        end
        count += 1
    end
    return jac_structures
end

function evaluateScalarQuadraticFunctionConstraintsJacobianStructure(constraints::Vector{OnePhase.ConstraintInfo{MathOptInterface.ScalarQuadraticFunction{Float64}}}, x::Array{Float64,1})
    jac_structures = nothing
    count = 1
    for constraint in constraints
        jac_structure = zeros(Float64, length(x))
        for term in constraint.func.affine_terms
            index = term.variable.value
            value = term.coefficient
            jac_structure[index] = value
        end
        if count == 1
            jac_structures = jac_structure
        else
            jac_structures = hcat(jac_structures, jac_structure)
        end
        count += 1
    end
    return jac_structures
end

function eval_a(m::Class_CUTEst, x::Array{Float64,1})
    if m.nlp != nothing
        a = cons(m.nlp, _cute_x(m, x) )
        if a == Float64[]
          a = zeros(1)
        end
        a_1 = a
        return [lb(a, m.bcon); ub(a, m.bcon); lb(x, m.bvar); ub(x, m.bvar)];
    end

    solver = m.solver
    # linear_constraints = vcat(solver.linear_le_constraints, solver.linear_ge_constraints, solver.linear_eq_constraints, solver.linear_int_constraints)
    # linear_constraints_evaluation = evaluateScalarAffineFunctionConstraints(linear_constraints, x)
    #
    # quadratic_constraints = vcat(solver.quadratic_le_constraints, solver.quadratic_ge_constraints, solver.quadratic_eq_constraints, solver.quadratic_int_constraints)
    # quadratic_constraints_evaluation = evaluateScalarQuadraticFunctionConstraints(quadratic_constraints, x)

    linear_constraints = vcat(solver.linear_eq_constraints, solver.linear_ge_constraints, solver.linear_le_constraints, solver.linear_int_constraints)
    linear_constraints_evaluation = evaluateScalarAffineFunctionConstraints(linear_constraints, x)

    quadratic_constraints = vcat(solver.quadratic_eq_constraints, solver.quadratic_ge_constraints, solver.quadratic_le_constraints, solver.quadratic_int_constraints)
    quadratic_constraints_evaluation = evaluateScalarQuadraticFunctionConstraints(quadratic_constraints, x)

    non_linear_constraints_evaluation = Float64[]
    if typeof(m.solver.nlp_data.evaluator) != EmptyNLPEvaluator
        non_linear_constraints_evaluation = zeros(Float64, length(m.solver.nlp_data.evaluator.constraints))
        MOI.eval_constraint(m.solver.nlp_data.evaluator, non_linear_constraints_evaluation, x)
    end

    a = vcat(linear_constraints_evaluation, quadratic_constraints_evaluation, non_linear_constraints_evaluation)
    if a == Float64[]
      a = zeros(1)
    end

    return [lb(a, m.bcon); ub(a, m.bcon); lb(x, m.bvar); ub(x, m.bvar)];
end

function _cute_x(m::Class_CUTEst, x::Array{Float64,1})
    if m.nlp != nothing
        ind = _i_not_fixed(m.nlp)
        #@show length(ind)
        if(length(x) != length(ind))
          error("$(length(x)) = length(x) != length(i_not_fixed) = $(length(ind))")
        end

        cute_x = deepcopy(m.nlp.meta.lvar) # get correct values of fixed variables
        cute_x[ind] = x

        return cute_x
    end

    ind = _i_not_fixed(m.solver.variable_info)
    #@show length(ind)
    if(length(x) != length(ind))
      error("$(length(x)) = length(x) != length(i_not_fixed) = $(length(ind))")
    end

    indexes = [i for i in 1:length(m.solver.variable_info)]
    lvar, uvar = extract_lvar_uvar(m.solver.variable_info, indexes)
    cute_x = deepcopy(lvar) # get correct values of fixed variables
    cute_x[ind] = x

    return cute_x
end

function eval_jac(m::Class_CUTEst, x::Array{Float64,1})
    if m.nlp != nothing
        cute_x = _cute_x(m, x)
        J_full_T = jac(m.nlp, cute_x)'
        J_T = J_full_T[_i_not_fixed(m.nlp),:];
        #my_eye = SparseArrays.speye(length(x))
        my_eye = SparseMatrixCSC{Float64}(LinearAlgebra.I, length(x), length(x))
        Q_T = [J_T[:,m.bcon.l_i] -J_T[:,m.bcon.u_i] my_eye[:,m.bvar.l_i] -my_eye[:,m.bvar.u_i]];
        return SparseArrays.sparse(Q_T')
    end

    cute_x = _cute_x(m, x)
    solver = m.solver
    linear_constraints = vcat(solver.linear_le_constraints, solver.linear_ge_constraints, solver.linear_eq_constraints, solver.linear_int_constraints)
    linear_constraints_jac_structure = evaluateScalarAffineFunctionConstraintsJacobianStructure(linear_constraints, x)

    quadratic_constraints = vcat(solver.quadratic_le_constraints, solver.quadratic_ge_constraints, solver.quadratic_eq_constraints, solver.quadratic_int_constraints)
    quadratic_constraints_jac_structure = evaluateScalarQuadraticFunctionConstraintsJacobianStructure(quadratic_constraints, x)
    d = m.solver.nlp_data.evaluator
    jac_struct = MOI.jacobian_structure(d)
    nnz_jac = length(jac_struct)
    J = zeros(Float64, nnz_jac)
    MOI.eval_constraint_jacobian(m.solver.nlp_data.evaluator, J, cute_x)

    J_full_T = nothing
    if linear_constraints_jac_structure!= nothing && !isempty(linear_constraints_jac_structure)
        J_full_T = hcat(linear_constraints_jac_structure)
        if quadratic_constraints_jac_structure != nothing && !isempty(quadratic_constraints_jac_structure)
            J_full_T = hcat(J_full_T, quadratic_constraints_jac_structure)
        end
        if J != nothing && !isempty(J)
            J_full_T = hcat(J_full_T, J)
        end
    elseif quadratic_constraints_jac_structure != nothing && !isempty(quadratic_constraints_jac_structure)
        J_full_T = hcat(quadratic_constraints_jac_structure)
        if J != nothing && !isempty(J)
            J_full_T = hcat(J_full_T, J)
        end
    elseif J != nothing && !isempty(J)
        count_variables = length(m.solver.variable_info)
        count_nonlinear_constraints = number_constraints_non_linear(m.solver)
        J_full_T = zeros(count_nonlinear_constraints, count_variables)
        for i in 1:length(jac_struct)
            J_full_T[jac_struct[i][1], jac_struct[i][2]] = J[i]
        end
        J_full_T = J_full_T'
    end

    J_T = J_full_T[_i_not_fixed(m.solver.variable_info),:];
    my_eye = SparseMatrixCSC{Float64}(LinearAlgebra.I, length(x), length(x))
    Q_T = [J_T[:,m.bcon.l_i] -J_T[:,m.bcon.u_i] my_eye[:,m.bvar.l_i] -my_eye[:,m.bvar.u_i]];
    return SparseArrays.sparse(Q_T')
end

function eval_grad_f(m::Class_CUTEst, x::Array{Float64,1})
    if m.nlp != nothing
        return grad(m.nlp, _cute_x(m, x))[_i_not_fixed(m.nlp)]
    end

    g = zeros(length(x))
    MOI.eval_objective_gradient(m.solver.nlp_data.evaluator, g, _cute_x(m, x))
    return g[_i_not_fixed(m.solver.variable_info)]
end

function get_reducedcosts(m::Class_CUTEst, y::Array{Float64,1})
    if m.nlp != nothing
        rc = zeros(m.nlp.meta.nvar)
        st_l = ncons_orginal(m)
        st_u = length(m.bvar.l) + ncons_orginal(m)
        rc[m.bvar.l_i] += y[(st_l+1):st_u]
        rc[m.bvar.u_i] -= y[(st_u+1):(length(m.bvar.u) + st_u)]
        return rc
    end

    rc = zeros(length(m.solver.variable_info))
    st_l = ncons_orginal(m)
    st_u = length(m.bvar.l) + ncons_orginal(m)
    rc[m.bvar.l_i] += y[(st_l+1):st_u]
    rc[m.bvar.u_i] -= y[(st_u+1):(length(m.bvar.u) + st_u)]
    return rc
end

function y_cons_net(m::Class_CUTEst, y::Array{Float64,1})
    return;
end

function eval_schur_and_jtprod(m::Class_CUTEst, x::Array{Float64,1}, y::Array{Float64,1}, s::Array{Float64,1})

end

function eval_lag_hess(m::Class_CUTEst, x::Array{Float64,1}, y::Array{Float64,1}, w::Float64)
    if m.nlp != nothing
        y_cons = zeros(m.nlp.meta.ncon)
        y_cons[m.bcon.l_i] -= y_l_con(y, m)
        y_cons[m.bcon.u_i] += y_u_con(y, m)
        H = hess(m.nlp, _cute_x(m, x), y_cons; obj_weight=w);

        H = LowerTriangular(H)
        ind = _i_not_fixed(m.nlp)
        H_not_fixed = H[ind,ind]
        return H_not_fixed
    end

    y_cons = zeros(number_constraints(m.solver))
    y_cons[m.bcon.l_i] -= y_l_con(y, m)
    y_cons[m.bcon.u_i] += y_u_con(y, m)

    d = m.solver.nlp_data.evaluator
    hessian_sparsity = MOI.hessian_lagrangian_structure(d)
    #Sparse hessian entry vector
    H_vec = zeros(Float64, length(hessian_sparsity))
    MOI.eval_hessian_lagrangian(d, H_vec, _cute_x(m, x), w, y_cons)
    H = spzeros(length(m.solver.variable_info), length(m.solver.variable_info))
    if !isempty(hessian_sparsity)
        index = 1
        if length(hessian_sparsity) < 3
            index = length(hessian_sparsity)
        elseif length(m.solver.variable_info) > 1
            index = 3
            for i in 3:length(m.solver.variable_info)
                index = index + i
            end
        end
        for i in 1:index
            H[hessian_sparsity[i][1], hessian_sparsity[i][2]]= H_vec[i]
        end
    end
    ind = _i_not_fixed(m.solver.variable_info)
    H_not_fixed = H[ind,ind]
    return H_not_fixed
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
function get_mutes_users_ids(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/mutes/users/ids.json", options))
end

function get_mutes_users_list(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/mutes/users/list.json", options))
end

ext: jl

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
__precompile__()

module Interact

using Reexport

@reexport using InteractBase
import InteractBase: notifications
import InteractBulma
import Widgets: Widget, @layout, @nodeps
import Observables: @on, @map!, @map

@reexport using DataStructures
@reexport using Observables
@reexport using Knockout
@reexport using CSSUtil
@reexport using WebIO
@reexport using Widgets

const notebookdir = joinpath(@__DIR__, "..", "doc", "notebooks")

const themes = Dict(
    :nativehtml => InteractBase.NativeHTML(),
    :bulma => InteractBulma.Bulma()
)

function InteractBase.settheme!(s::Symbol)
    (s in keys(themes)) || error("Theme $s is not supported")
    settheme!(themes[s])
end

function __init__()
    empty!(InteractBase.backend)
    settheme!(:bulma)
    nothing
end

end


text:
knn = KNNRegressor()
forest = RandomForestRegressor()
nusvr = NuSVR()

grid = Grid(resolution = 100, shuffle = false)

function fitknn(traindf; klower::Int = 1, kupper::Int = 20)
    folds = getfolds(traindf)
    r1 = range(knn, :K, lower = klower, upper = kupper)

    knnm = TunedModel(
        model = knn,
        resampling = folds,
        tuning = grid,
        range = [r1],
        measure = l2,
    )

    doublefit(knnm, traindf)
end

function fitnusvr(traindf, nuvalues, gammavalues)
    folds = getfolds(traindf)
    r_gamma = range(nusvr, :gamma; values = gammavalues)
    r_nu = range(nusvr, :nu, values = nuvalues)

    nusvrm = TunedModel(
        model = nusvr,
        resampling = folds,
        tuning = grid,
        range = [r_gamma, r_nu],
        measure = l2,
    )

    doublefit(nusvrm, traindf)
end

function fitforest(
    traindf;
    ntrees = [100, 200, 500],
    maxdepths = [-1, 10, 100],
    minleafs = [1, 5, 20, 100]
)
    folds = getfolds(traindf)
    r_n_trees =
        range(forest, :n_trees; values = ntrees)
    r_max_depth = range(forest, :max_depth, values = maxdepths)
    r_min_leaf = range(forest, :min_samples_leaf, values = minleafs)

    forestm = TunedModel(
        model = forest,
        resampling = CV(nfolds=5),
        tuning = grid,
        range = [r_n_trees, r_max_depth, r_min_leaf],
        measure = l2,
        acceleration = CPUProcesses(),
    )

    doublefit(forestm, traindf)
end

function fitcustommodel(traindf, model)
    folds = getfolds(traindf)
    doublefit(model, traindf)
end

function getfolds(data)
    fold = Tuple{Array{Int64,1},Array{Int64,1}}[]
    samples_tags = [1, 2, 3, 4, 5]
    for test_tag in samples_tags
        is_test = data.path_sample .== test_tag
        test = findall(is_test)
        train = findall(.!(is_test))
        push!(fold, (train, test))
    end

    return fold
end

function doublefit(model, traindf)
    X, y = preparedata(traindf)
    m_lon = machine(model, X, y[:, 1])
    m_lat = machine(model, X, y[:, 2])
    fit!(m_lon)
    fit!(m_lat)

    return m_lon, m_lat
end

function preparedata(data)
    X, y = unpack(
        data,
        X -> X ∈ [:magnetometr_x, :magnetometr_y, :magnetometr_z],
        y -> y ∈ [:lon, :lat],
    )
    X = coerce(
        X,
        :magnetometr_x => Continuous,
        :magnetometr_y => Continuous,
        :magnetometr_z => Continuous,
    )
    y = coerce(y, :lon => Continuous, :lat => Continuous)

    return X, y
end

function testmachines(machine_lon, machine_lat, testdf, coordsstand)

    test_X, test_y = preparedata(testdf)
    test_y[!, [:lon, :lat]] = MLJ.inverse_transform(coordsstand, test_y[!, [:lon, :lat]])

    ȳlon = predict(machine_lon, test_X)
    ȳlat = predict(machine_lat, test_X)
    predy = DataFrame(lon=ȳlon, lat=ȳlat)
    predy = coerce(predy, :lon => Continuous, :lat => Continuous)
    predy = MLJ.inverse_transform(coordsstand, predy)

    ȳperdicted = collect(zip(predy[:, :lon], predy[:, :lat]))
    ytest = collect(zip(test_y[:, :lon], test_y[:, :lat]))

    evaluateresults(ytest, ȳperdicted)
end

function evaluateresults(ytest, ypred)
    error = vechaversine(ytest, ypred)

    return Dict(
        "mean" => mean(error),
        "max" => maximum(error),
        "min" => minimum(error),
        "median" => median(error),
        "std" => std(error),
        "80" => quantile(error, 0.8),
    )
end

function vechaversine(x, y)
    if length(x) != length(y)
        throw(ArgumentError("x and y lengths must be equal"))
    end

    r = 6371000
    result = []
    for i = 1:length(x)
        point1 = x[i]
        point2 = y[i]
        push!(result, haversine(point1, point2, r))
    end

    return result
end

max_stars_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_issues_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_forks_repo_name: SebastianSuchowiak/MagneticLocSuchowiak


text:
# x, y without type, this is a generic function
function f(x, y)
    x + y
end

println(f(1, 2))
# 3
g = f # g is f
println(g(1, 2))

function with_return()
    return "I have return"
end
println(with_return())

function hypot(x, y)
    x = abs(x)
    y = abs(y)
    if x > y
        r = y/x
        return x*sqrt(1+r*r)
    end
    if y == 0
        return zero(x)
    end
    r = x/y
    return y*sqrt(1+r*r)
end

println(hypot(3, 4))
# 5

println(+(1, 2, 3) == (1 + 2 + 3))
# operator is function

max_stars_repo_name: dannypsnl/languages-learn

max_issues_repo_name: dannypsnl/languages-learn

max_forks_repo_name: dannypsnl/languages-learn


text:
module SpringCollab2020CustomBirdTutorial

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/CustomBirdTutorial" CustomBirdTutorial(x::Integer, y::Integer,
    birdId::String="birdId", onlyOnce::Bool=false, caw::Bool=true, faceLeft::Bool=true, info::String="TUTORIAL_DREAMJUMP", controls::String="DownRight,+,Dash,tinyarrow,Jump")

const placements = Ahorn.PlacementDict(
    "Custom Bird Tutorial (Spring Collab 2020)" => Ahorn.EntityPlacement(
        CustomBirdTutorial
    )
)

Ahorn.editingOptions(entity::CustomBirdTutorial) = Dict{String, Any}(
    "info" => String["TUTORIAL_CLIMB", "TUTORIAL_HOLD", "TUTORIAL_DASH", "TUTORIAL_DREAMJUMP", "TUTORIAL_CARRY", "hyperjump/tutorial00", "hyperjump/tutorial01"]
)

sprite = "characters/bird/crow00"

function Ahorn.selection(entity::CustomBirdTutorial)
    x, y = Ahorn.position(entity)
    scaleX = get(entity.data, "faceLeft", true) ? -1 : 1

    return Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y, sx=scaleX, jx=0.5, jy=1.0)]
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomBirdTutorial, room::Maple.Room)
    scaleX = get(entity.data, "faceLeft", true) ? -1 : 1

    Ahorn.drawSprite(ctx, sprite, 0, 0, sx=scaleX, jx=0.5, jy=1.0)
end

end


text:
"""
```
one_draw(m::AbstractDSGEModel, data::Matrix{Float64}; use_chand_recursion::Bool = true,
         verbose::Symbol = :low)
```
Finds and returns one valid draw from parameter distribution, along with its log likelihood and log posterior.
"""
function one_draw(m::AbstractDSGEModel, data::Matrix{Float64};
                  use_chand_recursion::Bool = true, verbose::Symbol = :low)
    success    = false
    draw       = vec(rand(m.parameters, 1))
    draw_loglh = draw_logpost = 0.0

    while !success
        try
            update!(m, draw)
            draw_loglh = likelihood(m, data, catch_errors = true,
                                    use_chand_recursion = use_chand_recursion,
                                    verbose = verbose)
            draw_logpost = prior(m)
            if (draw_loglh == -Inf) | (draw_loglh === NaN)
                draw_loglh = draw_logpost = -Inf
            end
        catch err
            if isa(err, ParamBoundsError)
                draw_loglh = draw_logpost = -Inf
            elseif isa(err, PosDefException) || isa(err, SingularException) ||
                   isa(err, LinearAlgebra.LAPACKException) || isa(err, DSGE.SteadyStateConvergenceError)
                draw_loglh = draw_logpost = -Inf
            else
                throw(err)
            end
        end
        if any(isinf.(draw_loglh))
            draw = vec(rand(m.parameters, 1))
        else
            success = true
        end
    end
    return vector_reshape(draw, draw_loglh, draw_logpost)
end

"""
```
initial_draw!(m::AbstractDSGEModel, data::Matrix{Float64}, c::ParticleCloud)
initial_draw!(m::AbstractDSGEModel, data::Matrix{Float64}, c::Cloud)
```

Draw from a general starting distribution (set by default to be from the prior) to
initialize the SMC algorithm. Returns a tuple (logpost, loglh) and modifies the
particle objects in the particle cloud in place.
"""
function initial_draw!(m::AbstractModel, data::Matrix{Float64},
                       c::Union{Cloud, ParticleCloud};
                       parallel::Bool = false, use_chand_recursion::Bool = true,
                       verbose::Symbol = :low)
    n_parts = length(c)

    # ================== Define closure on one_draw function ==================
    sendto(workers(), m = m)
    sendto(workers(), data = data)
    sendto(workers(), verbose = verbose)
    sendto(workers(), use_chand_recursion = use_chand_recursion)

    one_draw_closure() = one_draw(m, data; use_chand_recursion = use_chand_recursion,
                                  verbose = verbose)
    @everywhere one_draw_closure() = one_draw(m, data;
                                              use_chand_recursion = use_chand_recursion,
                                              verbose = verbose)
    # =========================================================================

    # For each particle, finds valid parameter draw and returns likelihood & posterior
    draws, loglh, logpost = if parallel
        @sync @distributed (vector_reduce) for i in 1:n_parts
            one_draw_closure()
        end
    else
        vector_reduce([one_draw_closure() for i in 1:n_parts]...)
    end

    update_draws!(c, draws)
    update_loglh!(c, vec(loglh))
    update_logpost!(c, vec(logpost))
    update_old_loglh!(c, zeros(n_parts))
    # Need to call `set_weights` as opposed to `update_weights`
    # since update_weights will multiply and 0*anything = 0
    set_weights!(c, ones(n_parts))
end

"""
```
function draw_likelihood(m, data, draw_vec; verbose::Symbol = :low)
```
Computes likelihood of a particular parameter draw; returns loglh and logpost.
"""
function draw_likelihood(m::AbstractDSGEModel, data::Matrix{Float64},
                         draw_vec::Vector{Float64}; verbose::Symbol = :low)
    update!(m, draw_vec)
    loglh   = likelihood(m, data, verbose = verbose)
    logpost = prior(m)
    return scalar_reshape(loglh, logpost)
end

"""
```
initialize_likelihoods!(m::AbstractDSGEModel, data::Matrix{Float64},
                        c::Union{Cloud, ParticleCloud};
                        parallel::Bool = false, verbose::Symbol = :low)
```
This function is made for transfering the log-likelihood values saved in the
Cloud from a previous estimation to each particle's respective old_loglh
field, and for evaluating/saving the likelihood and posterior at the new data, which
here is just the argument, data.
"""
function initialize_likelihoods!(m::AbstractDSGEModel, data::Matrix{Float64},
                                 c::Union{Cloud, ParticleCloud};
                                 parallel::Bool = false, verbose::Symbol = :low)
    n_parts = length(c)
    draws = (typeof(c) <: Cloud) ? get_vals(c; transpose = false) : Matrix{Float64}(get_vals(c)')

    # Retire log-likelihood values from the old estimation to the field old_loglh
    update_old_loglh!(c, get_loglh(c))

    # ============== Define closure on draw_likelihood function ===============
    sendto(workers(), m = m)
    sendto(workers(), data = data)
    sendto(workers(), verbose = verbose)

    draw_likelihood_closure(draw::Vector{Float64}) = draw_likelihood(m, data, draw;
                                                                 verbose = verbose)
    @everywhere draw_likelihood_closure(draw::Vector{Float64}) = draw_likelihood(m, data,
                                                                 draw; verbose = verbose)
    # =========================================================================

    # TODO: handle when the likelihood with new data cannot be evaluated (returns -Inf),
    # even if the likelihood was not -Inf prior to incorporating new data
    loglh, logpost = if parallel
        @sync @distributed (scalar_reduce) for i in 1:n_parts
            draw_likelihood_closure(draws[i, :])
        end
    else
        scalar_reduce([draw_likelihood_closure(draws[i, :]) for i in 1:n_parts]...)
    end
    update_loglh!(c, loglh)
    update_logpost!(c, logpost)
end

"""
```
function initialize_cloud_settings!(m::AbstractDSGEModel, cloud::ParticleCloud;
                                    tempered_update::Bool = false)
```
Initializes stage index, number of Φ stages, c, resamples, acceptance, and sampling time.
"""
function initialize_cloud_settings!(m::AbstractDSGEModel,
                                    cloud::Union{ParticleCloud,Cloud};
                                    tempered_update::Bool = false)
    if tempered_update
        cloud.ESS = [cloud.ESS[end]]
    else
        cloud.ESS[1] = get_setting(m, :n_particles)
    end
    cloud.stage_index = 1
    cloud.n_Φ         = get_setting(m, :n_Φ)
    cloud.resamples   = 0
    cloud.c           = get_setting(m, :step_size_smc)
    cloud.accept      = get_setting(m, :target_accept)
    cloud.total_sampling_time = 0.
    cloud.tempering_schedule  = zeros(1)
end


text:
# ------------ Solo map ------------ #

mutable struct SoloMap{K} <: AddressMap{K}
    pointer::AddressMap{<: K}
    SoloMap{K}() where K = new{K}()
    SoloMap(pointer::K) where K = new{K}(pointer)
end

@inline shallow_iterator(sm::SoloMap{K}) where K = NVector((:pointer, sm.pointer))

@inline function get(sm::SoloMap{Value}, addr, fallback)
    hasfield(sm, addr) || return fallback
    return getfield(sm, addr)
end

@inline function get_sub(sm::SoloMap, addr::A) where A <: Address
    hasfield(sm, addr) && return getfield(sm, addr)
    Empty()
end
@inline get_sub(sm::SoloMap, addr::Tuple{}) = Empty()

@inline getindex(sm::SoloMap, addrs...) = get_value(get_sub(sm, addrs))

@inline Base.isempty(sm::SoloMap) = isempty(getfield(sm, :pointer))

@inline haskey(sm::SoloMap, addr::A) where A <: Address = hasfield(sm, addr) && has_value(get_sub(sm, addr))

@inline has_sub(sm::SoloMap, addr::A) where A <: Address = hasfield(sm, addr)

function set_sub!(sm::SoloMap{K}, addr::A, v::AddressMap{<: K}) where {K, A <: Address}
    hasfield(sm, addr) || return
    !isempty(v) && setfield!(sm, addr, v)
end
function set_sub!(sm::SoloMap{K}, addr::A, v::V) where {K, V, A <: Address}
    hasfield(sm, addr) || return
    !isempty(v) && setfield!(sm, addr, convert(K, v))
end

function merge(sel1::SoloMap{K}, sel2::SoloMap{K}) where K
    new = SoloMap{K}()
    setfield!(new, :pointer, sel2.pointer)
    new, true
end
merge(sm::SoloMap, ::Empty) = SoloMap(sm.pointer), false
merge(::Empty, sm::SoloMap) = SoloMap(sm.pointer), false

function merge!(sel1::SoloMap{K}, sel2::SoloMap{K}) where K
    setfield!(sel1, :pointer, sel2.pointer)
    true
end
merge!(sm::SoloMap, ::Empty) = sm, false
merge!(::Empty, sm::SoloMap) = sm, false

function collect!(par::T, addrs::Vector{Any}, chd::Dict{Any, Any}, tr::S, meta) where {T <: Tuple, S <: SoloMap}
    iterate(tr) do (k, v)
        collect!((par..., k), addrs, chd, v, meta)
    end
end
function collect!(addrs::Vector{Any}, chd::Dict{Any, Any}, tr::S, meta) where S <: SoloMap
    iterate(tr) do (k, v)
        collect!((k, ), addrs, chd, v, meta)
    end
end

# Select.
function select(sm::SoloMap{K}) where K
    new = DynamicTarget()
    for (k, v) in shallow_iterator(sm)
        set_sub!(new, k, select(v))
    end
    new
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "rr"
version = v"5.5"

# Collection of sources required to build rr
sources = [
    GitSource("https://github.com/Keno/rr.git",
              "43288faa052adeed277e32b22611b3e649437e3b")
]

# Bash recipe for building across all platforms
script = raw"""
pip3 install pexpect
cd ${WORKSPACE}/srcdir/rr/

mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release \
      -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -Ddisable32bit=ON -DBUILD_TESTS=OFF -DWILL_RUN_TESTS=OFF -Dstaticlibs=ON ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# rr only supports Linux
platforms = [
    Platform("x86_64", "linux", libc="glibc"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("rr", :rr),
]

# Dependencies that must be installed before this package can be built
# This is really a build dependency
dependencies = [
    BuildDependency("capnproto_jll"),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
               preferred_gcc_version=v"8")

max_stars_repo_name: AbhimanyuAryan/Yggdrasil

max_issues_repo_name: AbhimanyuAryan/Yggdrasil

max_forks_repo_name: AbhimanyuAryan/Yggdrasil


text:
module TestSpecDef

export TestSpec, DataType, A, AL, ALB, DAA, DAAL, DAALS, DAALB, DAAS20, D, DS, DD, ND

@enum DataType A AL ALB DAA DAAL DAALS DAALB DAAS20 D DS DD ND

struct  TestSpec
    dataTypeID::DataType
    problemNbr::Int
    benchMarkSeconds::Int
    benchMarkSamples::Int
    dataDict::String
    directMode::Bool
    optimization::Bool
    writeToLpFile::Bool
    writeToExcel::Bool
    excelFile::String
    excelLocation::String
    sheetNbr::Int
end

end

max_stars_repo_name: SonntagK/JuMPContainerBenchmarks

max_issues_repo_name: SonntagK/JuMPContainerBenchmarks

max_forks_repo_name: SonntagK/JuMPContainerBenchmarks


text:
module MultiFieldFEOperatorsTests

using Gridap.Arrays
using Gridap.Algebra
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.Fields
using Gridap.Integration
using Gridap.MultiField
using Gridap.CellData
using SparseArrays
using Test
using LinearAlgebra

using Gridap.MultiField

order = 2

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

trian = get_triangulation(model)
degree = order
quad = CellQuadrature(trian,degree)

strian = SkeletonTriangulation(model)
sdegree = order
squad = CellQuadrature(strian,sdegree)

V = TestFESpace(model=model,order=order,reffe=:Lagrangian,conformity=:H1,valuetype=Float64)
Q = TestFESpace(model=model,order=order-1,reffe=:Lagrangian,conformity=:L2,valuetype=Float64)

U = TrialFESpace(V)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

function a(x,y)
  u,p = x
  v,q = y
  v*u + v*p - q*p
end

function l(y)
  v,q = y
  v*4 + q
end

function a_Γ(x,y)
  u,p = x
  v,q = y
  jump(v)*mean(u) + jump(∇(q))⋅jump(∇(p)) - mean(v)*mean(p)
end

t_Ω = AffineFETerm(a,l,trian,quad)
t_Γ = LinearFETerm(a_Γ,strian,squad)

op = AffineFEOperator(X,Y,t_Ω,t_Γ)

op = AffineFEOperator(SparseMatrixCSR,X,Y,t_Ω,t_Γ)

op = FEOperator(X,Y,t_Ω,t_Γ)
xh = zero(X)
b = residual(op,xh)
A = jacobian(op,xh)
test_fe_operator(op,get_free_values(xh),b)

function r(x,y)
  u,p = x
  v,q = y
  v*(u*u) + v*p*u - q*p - v*4 + q
end

function j(x,dx,y)
  u,p = x
  du,dp = dx
  v,q = y
  2*v*u*du + v*dp*u + v*p*du - q*dp
end

t_Ω = FETerm(r,j,trian,quad)
t_Ω_auto = FETerm(r,trian,quad)

x = FEFunction(X,rand(num_free_dofs(X)))
dx = get_cell_basis(X)
y = get_cell_basis(Y)

cell_r = get_cell_residual(t_Ω,x,y)
cell_j = get_cell_jacobian(t_Ω,x,dx,y)

cell_r_auto = get_cell_residual(t_Ω_auto,x,y)
cell_j_auto = get_cell_jacobian(t_Ω_auto,x,dx,y)
test_array(cell_r_auto,cell_r,≈)
test_array(cell_j_auto,cell_j,≈)

op = FEOperator(X,Y,t_Ω_auto)
xh = zero(X)
b = residual(op,xh)
A = jacobian(op,xh)
test_fe_operator(op,get_free_values(xh),b)

end # module


text:
using DashBootstrapComponents, DashHtmlComponents

number_input = html_div(
    [
        html_p("Type a number outside the range 0-10"),
        dbc_input(type="number", min=0, max=10, step=1),
    ],
    id="styled-numeric-input",
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
"""
    mstr(m::MarkerStyle)

Internal function to help in the specific case where a line with markers of
different color than the line is required.
See `add_sub_marker!` below.
"""
mstr(m::MarkerStyle) = "mk_$(m.marker)_$(col2str(m.color; str=true))"

"""
    add_sub_marker!(f, m)

Internal function to add a subroutine to the GLE script to deal with markers
that  must have a different color than the line they are associated with.
For instance if you want a blue line with red markers, you need to define a
specific subroutine for red-markers otherwise both line and markers are going
to be of the same color.
"""
function add_sub_marker!(
            f::Figure,
            m::MarkerStyle
        )::Nothing

    mstr(m) ∈ keys(f.subroutines) && return
    f.subroutines[str(m)] = """
        sub _$(mstr(m)) size mdata
        gsave
        set color $(col2str(m.color))
        marker $(m.marker) size
        grestore
        end sub
        define marker $(mstr(m)) _$(mstr(m))
        """
    return
end


text:
using Base.Filesystem:ispath
using Printf:@sprintf
using Plots
function download_dataset(filename)
    url = "http://elib.zib.de/pub/mp-testdata/tsp/tsplib/tsp/"*filename
    #println(url)
    if !isfile(filename)
        download(url,filename)
    else
        message = @sprintf("The file %s is exists!",filename)
        @info message
    end
end
function get_data(filename)
    data_list = []
    open(filename,"r") do stream
        flag = false
        for line in eachline(stream)
            if occursin("EOF",line)
                break
            elseif occursin("NODE_COORD_SECTION",line)
                flag = true
            elseif flag
                tp_list = Int64[]
                output = split(line,['\t',' '])
                vala = parse(Float32,output[1])
                valb = parse(Float32,output[2])
                valc = parse(Float32,output[3])
                append!(tp_list,vala)
                append!(tp_list,valb)
                append!(tp_list,valc)
                append!(data_list,tp_list)
            end
        end
        return data_list
    end
    arr_len = Int64(length(data_list)/3)
    data_list = reshape(data_list,3,arr_len)
    return data_list'
end
function main()
    filename = "eil51.tsp"
    download_dataset(filename)
    citys_mat = get_data(filename)
    # 初始化参数
    t0, tf, alpha, markov_len = 100, 1, 0.99, 10000
    # 获取城市的个数
    citys_num = size(citys_mat)[1]
    # 获取邻接矩阵
    citys_x = reshape(citys_mat[:,2],citys_num,1).*ones(1,citys_num)
    citys_y = reshape(citys_mat[:,3],citys_num,1).*ones(1,citys_num)
    citys_index = citys_mat[:,1]
    citys_distance = sqrt.((citys_x-citys_x').^2+(citys_y-citys_y').^2)
    E_current = Inf32 # 当前解对应的距离
    E_best = Inf32 # 表示的是最优解
    sol_new = [i for i in range(1,stop=citys_num,step=1)] # 初始化解
    sol_current =copy(sol_new)
    sol_best = copy(sol_new)
    len_list = Float32[]
    t = t0
    epoches = 0
    while t>=tf
        for step=1:markov_len
            # 产生随机扰动
            if rand() < 0.5
                # 两交换
                ind1,ind2 = 0,0
                while (ind1 == ind2)
                    ind1 = Int(ceil(rand()*citys_num))
                    ind2 = Int(ceil(rand()*citys_num))
                    sol_new[ind1],sol_new[ind2] = sol_new[ind2],sol_new[ind1]
                end
            else
                ind1,ind2,ind3 = 0,0,0
                while (ind1 == ind2) || (ind1==ind3)|| (ind2==ind3)
                    ind1 = Int(ceil(rand()*citys_num))
                    ind2 = Int(ceil(rand()*citys_num))
                    ind3 = Int(ceil(rand()*citys_num))
                end
                tp_list = [ind1,ind2,ind3]
                tp_list = sort(tp_list)
                ind1,ind2,ind3 = tp_list[1],tp_list[2],tp_list[3]
                tpa = copy(sol_new[ind1:ind2-1])
                tpb = copy(sol_new[ind2:ind3-1])
                sol_new[ind1:ind1+ind3-ind2-1],sol_new[ind1+ind3-ind2:ind3-1] = tpb,tpa
            end
            # 计算目标函数值
            E_new = 0
            for k in range(1,stop=citys_num-1,step=1)
                E_new+= citys_distance[sol_new[k],sol_new[k+1]]
            end
            E_new += citys_distance[sol_new[1],sol_new[end]]
            if E_new <E_current
                E_current = E_new
                sol_current = copy(sol_new)
                if E_new < E_best
                    E_best = E_new
                    sol_best = copy(sol_new)
                end
            else
                if rand()<exp(-(E_new-E_current)/t)
                    E_current = E_new
                    sol_current = copy(sol_new)
                else
                    sol_new = copy(sol_current)
                end
            end
        end
        append!(len_list,E_best)
        epoches += 1
        message = @sprintf("epoches: %d,temperature:%.2f,length:%0.4f",epoches,t,E_best)
        if E_best<0.1
            break
        end
        t = t*alpha
        @info message
    end
    #draw(len_list,sol_best,citys_mat)
end
#=
function draw(len_list,points_list,citys_mat)
    #plot(len_list,linewidth=2,title="The best length")
    #savefig("length.png")
    #closeall()
    len_arr = length(points_list)
    for k in range(1,stop=len_arr-1,step=1)
        start_id = points_list[k]
        end_id = points_list[k+1]
        tp_x = [citys_mat[start_id,2],citys_mat[end_id,2]]
        tp_y = [citys_mat[start_id,3],citys_mat[end_id,3]]
        plot!(tp_x,tp_y,linewidth=2)
    end
    start_id = points_list[end]
    end_id = points_list[1]
    tp_x = [citys_mat[start_id,2],citys_mat[end_id,2]]
    tp_y = [citys_mat[start_id,3],citys_mat[end_id,3]]
    title!("The path of the route")
    savefig("points.png")
    closeall()
end
=#
main()

max_stars_repo_name: MobTgZhang/CIAlgorithms

max_issues_repo_name: MobTgZhang/CIAlgorithms

max_forks_repo_name: MobTgZhang/CIAlgorithms


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

using Test
using TimeSeriesEcon

include("test_mit.jl")
include("test_tseries.jl")
include("test_mvtseries.jl")
include("test_workspace.jl")


text:
using VisualRegressionTests
using Plots

@testset "utils" begin
    gr()

    # ##########
    # # const. #
    # ##########

    max_x = 3e-6
    max_y = 10e-6
    max_t = 0.1e-12

    nx = 60
    ny = 200

    λ = 2.04e-6

    ϵ = 9.
    μ = 1.

    ϵ_defect = 1.
    xs_defect = [0, 1e-6, -1e-6]
    ys_defect = [1e-6, 2e-6, 3e-6]
    rs_defect = [0.5e-6, 0.1e-6, 0.2e-6]

    # ##############
    # # components #
    # ##############

    grid = Grid(max_x, max_y, max_t, nx, ny)
    light = Light(λ)
    permittivity = Permittivity(ϵ, grid)
    permeability = Permeability(μ, grid)

    implant!(permittivity, ϵ_defect, xs_defect, ys_defect, rs_defect, grid)

    # #############
    # # simulator #
    # #############
    s = Simulator(grid, light, permittivity, permeability)
    simulate!(s)

    @plottest begin
        plot_ϵ(s)
    end "assets/permittivity.png"

    @plottest begin
        plot_e_field(s)
    end "assets/e_field.png"
end


text:
using Test
using TrackTools
using TrackTools: interpolate_position
using TrackTools: matrix_line, calc_intersect_coordinates, hasintersect
using TrackTools: interpolate_time, find_intersections, TrackPoint
using TimeZones
using KMLTracks

@testset "interpolate_position" begin
    x1, y1, t1 = 0.0, 0.0, ZonedDateTime(2019, 11, 1, 10, 0, 0, tz"UTC")
    x2, y2, t2 = 2.0, 1.0, ZonedDateTime(2019, 11, 1, 10, 0, 40, tz"UTC")
    @test interpolate_position(x1, y1, t1, x2, y2, t2, ZonedDateTime(2019, 11, 1, 10, 0, 20, tz"UTC")) == (1.0, 0.5)
end

@testset "intersect sample" begin
    # positions
    x1, y1, t1 = 0.0, 0.2, ZonedDateTime(2019, 11, 1, 10, 0, 0, tz"UTC")
    x2, y2, t2 = 1.0, 0.2, ZonedDateTime(2019, 11, 1, 10, 0, 40, tz"UTC")

    # line
    x1L, y1L = 0.5, -0.2
    x2L, y2L = 0.5, 0.4

    # calculate matrix for line
    AL = matrix_line(x1L, y1L, x2L, y2L)

    # calculate coordinates of intersect
    xI, yI = calc_intersect_coordinates(x1, y1, x2, y2, AL)

    @test (xI, yI) == (0.5, 0.2)

    @test hasintersect(x1, y1, x2, y2, xI, yI)

    t = interpolate_time(x1, y1, t1, x2, y2, t2, xI, yI)
    @test t == ZonedDateTime(2019, 11, 1, 10, 0, 20, tz"UTC")

end

@testset "intersect sample" begin
    fname = "sample.kml"
    kmldoc = read_kml_file(fname)
    @test length(kmldoc.placemark.track.points) == 45
    #y1L, x1L = (46.585904, 0.309138)  # LBFI ALPHA2 SUD (lat/lon)
    #y2L, x2L = (46.585963, 0.309425)  # LFBI ALPHA2 NORD (lat/lon)

    y1L, x1L = (46.586296, 0.308861)  # LFBI ALPHA SUD
    y2L, x2L = (46.586480, 0.309045)  # LFBI ALPHA NORD

    points = ((p.time, p.longitude, p.latitude) for p in kmldoc.placemark.track.points)
    lst_I = find_intersections(x1L, y1L, x2L, y2L, points)

    @test length(lst_I) == 1
    @test lst_I[1].time == ZonedDateTime(2019, 11, 7, 8, 44, 49, 473, tz"UTC")

end


text:
#used
function storage_costs(data,year::Int64=2021)
    if (year==2021)
        #csts=[944,472,236,118,59,29.5,14.75,7.375]
        #csts=[472,236,118,59,29.5,14.75,7.375]
        #csts=[236,118,59,29.5,14.75,7.375]
        csts=[29.5,14.75,7.375,3.6875,1.84375]
    elseif (year==2030)
        #csts=[560,280,140,70,35,17.5,8.75,4.375]
        #csts=[280,140,70,35,17.5,8.75,4.375]
        x=1
        csts=[35/x,17.5/x,8.75/x,4.375/x,2.1875/x,1.09375/x]
        #csts=[1,1,1,1,1,1]
    else
        println("No battery cost data for specified year defaulting to 2021.")
        #csts=[944,472,236,118,59,29.5,14.75,7.375]
        #csts=[472,236,118,59,29.5,14.75,7.375]
        csts=[236,118,59,29.5,14.75,7.375]
    end
    [data["ne_storage"][string(i)]["inst_cost"]=data["ne_storage"][string(i)]["eq_cost"]=cst for (i,cst) in enumerate(csts)]
    [data["ne_storage"][string(i)]["eq_cost"]=cst for (i,cst) in enumerate(csts)]
end


text:
using Test
using ScottishTaxBenefitModel
using .ModelHousehold: Household, Person, People_Dict, is_single,
    default_bu_allocation, get_benefit_units, get_head, get_spouse, search,
    pers_is_disabled, pers_is_carer, printpids
using .ExampleHouseholdGetter
using .Definitions
using .Results: HousingResult
using .FRSHouseholdGetter
using .GeneralTaxComponents: WEEKS_PER_YEAR
using .RunSettings: Settings
using .LocalLevelCalculations: apply_size_criteria, apply_rent_restrictions,
    make_la_to_brma_map, LA_BRMA_MAP, lookup, apply_rent_restrictions, calc_council_tax,
    LA_NAMES, LA_CODES

using .STBParameters
using .Intermediate: make_intermediate, MTIntermediate
using .ExampleHelpers

## FIXME don't need both
lmt = LegacyMeansTestedBenefitSystem{Float64}()
sys = get_system( scotland=true )

rc = @timed begin
    num_households,total_num_people,nhh2 = FRSHouseholdGetter.initialise( DEFAULT_SETTINGS )
end

@testset "LHA and assoc. mappings" begin
    # basic test/retrieve
    # println( LA_BRMA_MAP )
    @test LA_BRMA_MAP.map[:S12000049] == :S33000009
    lmt
    @test lookup( sys.hr.brmas, :S12000049, 4 ) == 322.19
end

@testset "Rooms Restrictions" begin

    hh = get_example( cpl_w_2_children_hh )
    hh.bedrooms = 12 # set to a big number
    delete_child!( hh ) # start with 1

    println( hh.tenure )
    hh.tenure = Private_Rented_Unfurnished

    println( sys.hr.maximum_rooms )
    bus = get_benefit_units(hh)
    bu = bus[1]

    # single_parent_hh single_hh childless_couple_hh
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    println( "got nbeds as $nbeds " )
    oldnbeds = 0
    @test nbeds == 2 # so 1 bed for adults + 1 shared

    np = add_child!( hh, 11, Female )
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    println( "got nbeds as $nbeds " )
    oldnbeds = 0
    age = 4
    # base case: 2 children aged 2 and 5: different genders (sexes?)
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    @test nbeds == 2 # so 1 bed for adults + 1 shared

    sys.hr.maximum_rooms = 5 # add 1 so we can test a bit more`
    np = add_child!( hh, 11, Female )
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    @test nbeds == 3 # so 1 bed for adults + 1 shared + 1 for 11 yo

    np = add_child!( hh, 11, Male )
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    @test nbeds == 4 # so 1 bed for adults + 1 shared 11,2 yo male + 1 F

    np = add_child!( hh, 12, Male )
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    @test nbeds == 4 # so 1 bed for adults + 2 shared + 1 for 11 M and F

    np = add_child!( hh, 13, Female )
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    @test nbeds == 4 # so 1 bed for adults + 2 shared + 1 for 11 M and F

    np = add_child!( hh, 15, Female )
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    @test nbeds == 5 # same as above - max should kick in
    hh = get_example( cpl_w_2_children_hh )

    for i in 1:0
        age += 1
        sex = iseven(i) ? Male : Female
        np = add_child!( hh, age, sex )
        oldnbeds = nbeds
        intermed = make_intermediate(
            hh,
            sys.hours_limits,
            sys.age_limits,
            sys.child_limits )
        nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
        nc = num_children( hh )
    end

    hh = make_hh() # all at defaults
    head = get_head( hh )
    head.age = 20
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    println( "beds for under 35s $nbeds ")
    @test nbeds == 0 # single room
    head.age = 40

    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    nbeds = apply_size_criteria( hh, intermed.hhint, sys.hr )
    println( "beds for over 35s $nbeds ")
    @test nbeds == 1 # single room + bed for over 35
end

@testset "Local Housing Allowance" begin
    @test sys.hr.rooms_rent_reduction ≈ [0.14, 0.25]
    sys.hr.maximum_rooms = 4 # set this back to actual
    for (name,hh) in get_all_examples()
        println( "on hhld $name")
        hh.tenure = Private_Rented_Furnished
        hh.gross_rent = 300.0
        intermed = make_intermediate(
            hh,
            sys.hours_limits ,
            sys.age_limits,
            sys.child_limits )
        rr = apply_rent_restrictions( hh, intermed.hhint, sys.hr )
        println( rr )
    end
    # this hhld is in Glasgow
    #
    for tenure in [Private_Rented_Furnished, Council_Rented]
        for adults in 1:2
            for kids in 0:5
                for age in [30,40,70]
                    if adults == 1
                        hh = make_hh( adults=adults, children=kids, age=age, tenure=tenure, rent=500.0 )
                    else
                        hh = make_hh( adults=adults, children=kids, age=age, spouse_age=age, tenure=tenure, rent=500.0 )
                    end
                    intermed = make_intermediate(
                        hh,
                        sys.hours_limits ,
                        sys.age_limits,
                        sys.child_limits )
                    rr = apply_rent_restrictions( hh, intermed.hhint, sys.hr )
                    if adults == 1
                        if age == 70 && tenure == Council_Rented
                            # no bedroom tax for hhls all ads over pension age socially renting
                            @test intermed.hhint.all_pension_age
                            println( "intermed.all_pension_age $(intermed.hhint.all_pension_age)")
                            @test rr.allowed_rooms == hh.bedrooms
                        elseif kids == 0
                            if age == 30
                                @test rr.allowed_rooms == 0
                            else
                                @test rr.allowed_rooms == 1
                            end
                        elseif kids == 1
                            @test rr.allowed_rooms == 2 # you & the child, regardless of age
                        elseif kids == 2
                            @test rr.allowed_rooms == 3
                        elseif 2 < kids > 4
                            @test rr.allowed_rooms ∈ [3,4] # depends on ages
                        elseif kids > 4
                            @test rr.allowed_rooms == 4
                        end

                    else # 2 adults
                        if age == 70 && tenure == Council_Rented
                            # no bedroom tax for hhls all ads over pension age socially renting
                            @test intermed.hhint.all_pension_age
                            println( "intermed.all_pension_age $(intermed.hhint.all_pension_age)")
                            @test rr.allowed_rooms == hh.bedrooms
                        elseif kids == 0
                            @test rr.allowed_rooms == 1
                        elseif kids == 1
                            @test rr.allowed_rooms == 2 # you & the child, regardless of age
                        elseif kids == 2
                            @test rr.allowed_rooms ∈ [2,3] # this depends on whether kids can share
                        elseif 2 < kids > 4
                            @test rr.allowed_rooms ∈ [3,4] # depends on kids can share
                        elseif kids > 4
                            @test rr.allowed_rooms == 4
                        end
                    end
                end
            end # kids
        end # adults
    end # tenure

    for tenure in [Private_Rented_Furnished, Council_Rented]
        for adults in 1:2
            for kids in 0:5
                for age in [30,40,70]
                    if adults == 1
                        hh = make_hh( adults=adults, children=kids, age=age, tenure=tenure, rent=500.0 )
                    else
                        hh = make_hh( adults=adults, children=kids, age=age, spouse_age=age, tenure=tenure, rent=500.0 )
                    end
                    intermed = make_intermediate(
                        hh,
                        sys.hours_limits ,
                        sys.age_limits,
                        sys.child_limits )
                    rr = apply_rent_restrictions( hh, intermed.hhint, sys.hr )

                    if tenure == Private_Rented_Furnished
                        # GLASGOW 2020/1
                        allowed = [80.55 113.92 149.59 172.6 322.19]
                        @test rr.allowed_rent ≈ allowed[rr.allowed_rooms+1]
                        @test rr.allowed_rooms ∈ 0:4
                    else
                        if age != 70
                            @test rr.allowed_rooms ∈ 0:4
                        end
                        if rr.excess_rooms == 0
                            @test rr.allowed_rent == 500
                        elseif rr.excess_rooms == 1
                            @test rr.allowed_rent ≈ 500*(1-0.14)
                        else
                            @test rr.allowed_rent ≈ 500*(1-0.25)
                        end
                    end
                end
            end
        end
    end

    num_restricted = 0
    bedroom_tax = 0
    for hhno in 1:num_households
        hh = FRSHouseholdGetter.get_household( hhno )
        # TODO UPRATE
        if hhno % 500 == 0
            println( "on hhld $hhno")
        end
        intermed = make_intermediate(
            hh,
            sys.hours_limits ,
            sys.age_limits,
            sys.child_limits )
        rr = apply_rent_restrictions( hh, intermed.hhint, sys.hr )
        if rr.excess_rooms > 0
            num_restricted += hh.weight
            if is_social_renter( hh.tenure )
                bedroom_tax += hh.weight
            end
        end
    end
    println( "initial run: number with excess rooms $num_restricted bedroom tax $bedroom_tax" )
end

@testset "Council Tax" begin
    by_band = Dict{CT_Band, Real}()
    by_la = Dict()
    for c in instances( CT_Band )
        by_band[c] = 0.0
    end
    for c in LA_CODES
        by_la[c] = [0.0, 0.0]
    end
    println( by_la[:S12000019][1] )
    value = 0.0
    dwellings = 0.0
    for hhno in 1:num_households
        hh = FRSHouseholdGetter.get_household( hhno )
        intermed = make_intermediate(
            hh,
            sys.hours_limits ,
            sys.age_limits,
            sys.child_limits )
        println( "ct band $(hh.ct_band) council $(hh.council)")
        ct = calc_council_tax( hh, intermed.hhint, sys.loctax.ct )
        by_band[hh.ct_band] += hh.weight
        by_la[hh.council][2] += ct*hh.weight*WEEKS_PER_YEAR
        by_la[hh.council][1] += hh.weight
        value += ct*hh.weight*WEEKS_PER_YEAR
        dwellings +=hh.weight
    end
    println( "band,num dwellings")
    for c in instances( CT_Band )
        println( "$c,$(trunc(by_band[c]))")
    end

    println( "name,ccode,hhlds,raised,av")
    for c in LA_CODES
        name = LA_NAMES[c]
        raised = trunc( by_la[c][2] )
        hhlds = trunc( by_la[c][1] )
        av = trunc( by_la[c][2]/by_la[c][1])
        println( "$(name),$c,$hhlds,$raised,$av")
    end

    println( "total raised $(trunc(value/1_000_000))m pa before rebates")
    println( "dwellings $(trunc(dwellings)) ")
    println( "av per dwelling, before ctrebate $(trunc(value/dwellings))")
    hh = make_hh(adults=2)
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    println( "ct band $(hh.ct_band) council $(hh.council)")
    ct = calc_council_tax( hh, intermed.hhint, sys.loctax.ct )
    @test hh.ct_band == Band_B
    @test ct ≈ 1_078.00/WEEKS_PER_YEAR # glasgow 2020/1 CT band b per week
end


text:
function mlogit_with_Z(θ, X, Z, y)
    α = θ[1:end-1]
    γ = θ[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigα = [reshape(α,K,J-1) zeros(K)]

    P = exp.(X*bigα .+ Z*γ) ./ sum.(eachrow(exp.(X*bigα .+ Z*γ)))

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

    return loglike
end


text:
using BenchmarkTools
using Distributed

function calc_pi_by_gcd(N)
    s = 0
    for a in 1:N
        for b in 1:N
            s += ifelse(gcd(a,b)==1, 1, 0)
        end
        s
    end
    println(sqrt(6N^2/s))
end

#@time calc_pi_by_gcd(20000)

function calc_pi_by_gcd_parallel(N)
    """
    Note that @parallel has been renamed into @difrom Julia 1.0
    https://discourse.julialang.org/t/loaderror-parallel/15126
    """
    c = @distributed (+) for a in 1:N
        s = 0
        for b in 1:N
            s += ifelse(gcd(a,b)==1, 1, 0)
        end
        s
    end
    println(sqrt(6N^2/c))
end

addprocs()

#@time calc_pi_by_gcd_parallel(20000)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module problem119

description = """
Digit power sum

The number 512 is interesting because it is equal to the sum of its digits raised to some power: 5 + 1 + 2 = 8, and 83 = 512. Another example of a number with this property is 614656 = 284.
We shall define an to be the nth term of this sequence and insist that a number must contain at least two digits to have a sum.
You are given that a2 = 512 and a10 = 614656.

Find a30.
"""

# sum the digits of an integer by repeated div/mod by 10
function digitsum(n :: Integer)
  total = 0
  while n > 0
    total += n % 10
    n = div(n, 10)
  end
  return total
end

# Find all numbers with are powers of their digit sums, up to a maximum number of digits
# Insight: the digit sum for any number < 10^N can be at most N*9
#          so just consider powers of i from 2 to 9N, checking digitsums
function digitPowerSums(maxDigits = 15)
  results = BigInt[]
  xMax = BigInt(10)^maxDigits
  for i = 2 : 9 * maxDigits
    x = BigInt(i)
    while x < xMax
      if x > 10 && digitsum(x) == i
        push!(results, x)
      end
      x *= i
    end
  end
  sort!(results)
  return results
end

using Base.Test

@test digitsum(1234) == (1+2+3+4)
a = digitPowerSums()
@test a[2] == 512
@test a[10] == 614656

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using Test
using GGPlot2, DataFrames

df_string = DataFrame(x = randn(10), y=randn(10), z1=rand(["a"; "b"], 10))

plt = ggplot(df_string) + geom_point(aes(x=:x, y=:y, color=:z1))

@test GGPlot2.isggplot(plt)


text:
#=
Largely copied from NamedDims:
https://github.com/invenia/NamedDims.jl/blob/master/src/broadcasting.jl
=#

using Base.Broadcast:
    Broadcasted, BroadcastStyle, DefaultArrayStyle, AbstractArrayStyle, Unknown

struct KeyedStyle{S <: BroadcastStyle} <: AbstractArrayStyle{Any} end
KeyedStyle(::S) where {S} = KeyedStyle{S}()
KeyedStyle(::S, ::Val{N}) where {S,N} = KeyedStyle(S(Val(N)))
KeyedStyle(::Val{N}) where {N} = KeyedStyle{DefaultArrayStyle{N}}()

function KeyedStyle(a::BroadcastStyle, b::BroadcastStyle)
    inner_style = BroadcastStyle(a, b)
    inner_style isa Unknown ? Unknown() : KeyedStyle(inner_style)
end

Base.BroadcastStyle(::Type{<:KeyedArray{T,N,AT}}) where {T,N,AT} =
    KeyedStyle{typeof(BroadcastStyle(AT))}()
Base.BroadcastStyle(::KeyedStyle{A}, ::KeyedStyle{B}) where {A, B} = KeyedStyle(A(), B())
Base.BroadcastStyle(::KeyedStyle{A}, b::B) where {A, B} = KeyedStyle(A(), b)
Base.BroadcastStyle(a::A, ::KeyedStyle{B}) where {A, B} = KeyedStyle(a, B())
Base.BroadcastStyle(::KeyedStyle{A}, b::DefaultArrayStyle) where {A} = KeyedStyle(A(), b)
Base.BroadcastStyle(a::AbstractArrayStyle{M}, ::KeyedStyle{B}) where {B,M} = KeyedStyle(a, B())

using NamedDims: NamedDimsStyle
# this resolves in favour of KeyedArray(NamedDimsArray())
Base.BroadcastStyle(a::NamedDimsStyle, ::KeyedStyle{B}) where {B} = KeyedStyle(a, B())
Base.BroadcastStyle(::KeyedStyle{A}, b::NamedDimsStyle) where {A} = KeyedStyle(A(), b)

function unwrap_broadcasted(bc::Broadcasted{KeyedStyle{S}}) where {S}
    inner_args = map(unwrap_broadcasted, bc.args)
    Broadcasted{S}(bc.f, inner_args)
end
unwrap_broadcasted(x) = x
unwrap_broadcasted(x::KeyedArray) = parent(x)

function Broadcast.copy(bc::Broadcasted{KeyedStyle{S}}) where {S}
    inner_bc = unwrap_broadcasted(bc)
    data = copy(inner_bc)
    R = broadcasted_keys(bc)
    KeyedArray(data, map(copy, R))
end

function Base.copyto!(dest::AbstractArray, bc::Broadcasted{KeyedStyle{S}}) where {S}
    inner_bc = unwrap_broadcasted(bc)
    data = copyto!(keyless(dest), inner_bc)
    R = unify_keys(keys_or_axes(dest), broadcasted_keys(bc))
    KeyedArray(data, R)
end

broadcasted_keys(bc::Broadcasted) = broadcasted_keys(bc.args...)
function broadcasted_keys(a, bs...)
    a_r = broadcasted_keys(a)
    b_r = broadcasted_keys(bs...)
    unify_longest(a_r, b_r)
end
broadcasted_keys(a::AbstractArray) = keys_or_axes(a)
broadcasted_keys(a) = tuple()

#===== Unification, also used by map, hcat, etc. =====#

function unify_longest(short::Tuple, long::Tuple)
    length(short) > length(long) && return unify_longest(long, short)
    overlap = unify_keys(short, ntuple(i -> long[i], length(short)))
    extra = ntuple(i -> long[i + length(short)], length(long) - length(short))
    return (overlap..., extra...)
end
unify_longest(x::Tuple) = x
unify_longest(x,y,zs...) = unify_longest(unify_longest(x,y), zs...)

unify_keys(left::Tuple, right::Tuple) = map(unify_one, left, right)
unify_keys(left::Tuple) = left
unify_keys(left, right, more...) = unify_keys(unify_keys(left, right), more...)

unifiable_keys(left::Tuple, right::Tuple) = map(who_wins, left, right) isa Tuple{Vararg{AbstractArray}}

function unify_one(x::AbstractArray, y::AbstractArray)
    out = who_wins(x,y)
    out === nothing && throw(ArgumentError("key vectors must agree; got $x != $y"))
    out
end
unify_one(x,y,zs...) = unify_one(unify_one(x,y), zs...)

"""
    who_wins(axisranges(A,1), axisranges(B,1))
    who_wins(r, s, t, ...)

For broadcasting, but also `map` etc, this compares individual key vectors
and returns the one to keep.
In general they must agree `==`, and the simpler type will be returned
(e.g. `Vector + UnitRange -> UnitRange`).

However default key vectors `Base.OneTo(n)` are regarded as wildcards.
They need not agree with anyone, and are always discarded in favour of other types.

If the keys disagree it returns `nothing`.
Call `unify_one()` to have an error instead.
"""
who_wins(x,y,zs...) = who_wins(who_wins(x,y), zs...)

who_wins(r::AbstractVector, s::AbstractVector) = r === s ? r : r == s ? r : nothing
who_wins(vec::AbstractVector, ran::AbstractRange) = vec == ran ? ran : nothing
who_wins(ran::AbstractRange, vec::AbstractVector) = vec == ran ? ran : nothing
who_wins(r::AbstractRange, s::AbstractRange) = r == s ? r : nothing

who_wins(arr::AbstractVector, ot::Base.OneTo) = arr
who_wins(ot::Base.OneTo, arr::AbstractVector) = arr
who_wins(ot::Base.OneTo, ot′::Base.OneTo) = ot # else ambiguous
who_wins(arr::AbstractRange, ot::Base.OneTo) = arr # also to solve ambiguity
who_wins(ot::Base.OneTo, arr::AbstractRange) = arr


text:
### Main Program: parse args and call main() ###
using ArgParse

using WWURay

""" parse_cmdline
Parse the command line args to specify scene, camera, and image size """
function parse_cmdline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--scene", "-s"
            help="scene"
            arg_type=String
            default="teapot"
        "--camera", "-c"
            help="camera"
            arg_type=Int
            default=1
        "width"
            help="image width"
            arg_type=Int
        "height"
            help="image height"
            arg_type=Int
        "--samples"
            help="number of samples"
            arg_type=Int
        "--denoise"
            help="denoise resulting image"
            arg_type=Bool
    end

    return parse_args(s)
end

a = parse_cmdline()

WWURay.main(a["scene"], a["width"], a["height"], a["camera"], a["samples"], a["denoise"])

max_stars_repo_name: aedans/raytracer

max_issues_repo_name: aedans/raytracer

max_forks_repo_name: aedans/raytracer


text:
using FiniteHorizonPOMDPs

# run(`clear`)

@testset "finite_horizon" begin
    pp = simple_pp()
    me = PMDPs.PMDPe(pp)

    # Solve Value Iteration
    solver = ValueIterationSolver(max_iterations=1000, belres=1e-3, include_Q=true);
    VIPolicy = DiscreteValueIteration.solve(solver, me);

    fhsolver = FiniteHorizonSolver(false)

    # FH
    FHPolicy = POMDPs.solve(fhsolver, me)

    for s in states(me)
        # if !PMDPs.sale_impossible(me,s)
            @test FiniteHorizonPOMDPs.action(FHPolicy, s) == action(VIPolicy, s)
            # display("$(FiniteHorizonPOMDPs.stage_stateindex(me, s, 1)) - $s - $(PMDPs.product(me, s))")
        # end

        # if FiniteHorizonPOMDPs.action(FHPolicy, s) != action(VIPolicy, s)
        #     print(PMDPs.sale_impossible(me, s), " ")
        #     display(s)
        # end
    end

end


text:
include("load.jl")
FoldsPreviewTests.runtests()


text:
# Utility functions used throughout the package

using IterTools: partition  # to split iterators
import DataFrames: DataFrame
using uCSV  # for reading tabular data

# replacement for DataFrames readtable
function readtable(file)
    DataFrame(uCSV.read(file,
              trimwhitespace = true,
              typedetectrows = 10,
              header = 1,
              comment = '#'))
end

# Custom exceptions

"Error raised when a particular method needs to be implemented for some type"
struct NotImplementedError <: Exception end

# Duplicate-finding

""" Get indices of the unique items in `itr`
    For non-unique items, returns the highest index.

    ```jl
    julia> WaterData.unique_indices([1,1,4,7,4,4])
    3-element Array{Int64,1}:
     2
     6
     4
    ```
    """
unique_indices(itr) = indexin(collect(unique(itr)), collect(itr))

# Miscellaneous utility funcs

const VectorPair{T} = NTuple{2, Vector{T}}

""" Map a function `f` across rows of the array `A`

    ```jl
    julia> A = [1 2 3;
                4 5 6];
    julia> WaterData.maprows(sum, A)
    2x1 Array{Int64,2}:
      6
     15
    ```
    """
maprows(f, A::Matrix) = mapslices(f, A, 2)

""" Get an iterator of each adjacent pair in a sequence

    ```jl
    julia> seq = [1,2,3,4];
    julia> collect(WaterData.adjacentpairs(seq))
    3-element Array{Tuple{Int64,Int64},1}:
     (1,2)
     (2,3)
     (3,4)
    ```
    """
adjacentpairs(itr) = partition(itr, 2, 1)

"2D analogue of the cross product, *v* × *w* = v₁w₂ - v₂w₁"
crossprod2d(v::Vector, w::Vector) = v[1]*w[2] - v[2]*w[1]

"Is x outside the closed interval [a, b]?"
isoutside(x, a, b) = x < a || x > b

"Check if the line segments `AB` and `CD` cross"
function intersects(A::Vector, B::Vector, C::Vector, D::Vector)
    # AB is the line from p to p + r
    # CD is the line from q to q + s
    p = A
    r = B - A
    q = C
    s = D - C

    rxs = crossprod2d(r, s)
    qpxr = crossprod2d((q - p), r)
    t = crossprod2d((q - p), s) ./ rxs
    u = qpxr ./ rxs

    if (rxs ≠ 0) && (0 <= t <= 1) && (0 <= u <= 1)
        # cross or touch
        return true
    else
        # may be colinear or disjoint, but they don't cross
        return false
    end
end
intersects(AB::VectorPair, CD::VectorPair) = intersects(AB[1], AB[2], CD[1], CD[2])


text:
# linear operator implementation for Petsc matrix-free

"""
  Petsc matrix-free linear operator.

  **Public Fields**

   * none
"""
mutable struct PetscMatFreeLO <: AbstractPetscMatFreeLO
  A::PetscMat  # shell matrix
  xtmp::PetscVec  # these are shared with the PC if possible
  btmp::PetscVec
  ctx
  is_shared::Bool
  is_finalized::Bool

  nsolves::Int
  ntsolves::Int

  # MPI stuff
  comm::MPI.Comm
  myrank::Int
  commsize::Int
end

"""
  Outer constructor for PetscMatFreeLO

  **Inputs**

   * pc: a Petsc preconditioner of some kind
   * mesh
   * sbp
   * eqn
   * opts
"""
function PetscMatFreeLO(pc::Union{AbstractPetscMatPC, AbstractPetscMatFreePC},
                    mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  A = createPetscMatShell(mesh, sbp, eqn, opts)

  # set operations:
  fptr = cfunction(applyLinearOperator_wrapper, PetscErrorCode, (PetscMat, PetscVec, PetscVec))
  MatShellSetOperation(A, PETSc2.MATOP_MULT, fptr)

  fptr = cfunction(applyLinearOperatorTranspose_wrapper, PetscErrorCode, (PetscMat, PetscVec, PetscVec))
  MatShellSetOperation(A, PETSc2.MATOP_MULT_TRANSPOSE, fptr)

  pc2 = getBasePC(pc)
  if typeof(pc2) <: PetscMatPC  # if pc2 has Petsc vectors inside it
    xtmp = pc2.xtmp
    btmp = pc2.btmp
  else
    xtmp = createPetscVec(mesh, sbp, eqn, opts)
    btmp = createPetscVec(mesh, sbp, eqn, opts)
  end

  ctx = C_NULL  # this gets set later
  is_shared = false

  is_finalized = false
  nsolves = 0
  ntsolves = 0
  comm = eqn.comm
  myrank = eqn.myrank
  commsize = eqn.commsize

  # update the PC so it knows if the LO used its matrix or not
  getBasePC(pc).is_shared = is_shared

  return PetscMatFreeLO(A, xtmp, btmp, ctx, is_shared, is_finalized, nsolves,
                        ntsolves, comm, myrank, commsize)
end

function free(lo::PetscMatFreeLO)

  if !lo.is_finalized
    if lo.A.pobj != C_NULL
      PetscDestroy(lo.A)
      lo.A.pobj = C_NULL
    end

    if lo.xtmp.pobj != C_NULL
      PetscDestroy(lo.xtmp)
      lo.xtmp.pobj = C_NULL
    end

    if lo.btmp.pobj != C_NULL
      PetscDestroy(lo.btmp)
      lo.btmp.pobj = C_NULL
    end
  end

  lo.is_finalized = true

  return nothing
end

"""
  This function sets the ctx for the underlying Petsc matrix object.  The
  end result is that the mesh, sbp, eqn, and opts that are passed into this
  function will be passed to the user-defined [`applyLinearOperator`](@ref)
  during the next linear solve.  See that function for a description of
  the arguments.

  **Inputs**

   * lo: the user defined linear operator
   * mesh
   * sbp
   * eqn
   * opts
   * ctx_residual
   * t
"""
function setLOCtx(lo::AbstractPetscMatFreeLO, mesh::AbstractMesh,
                  sbp::AbstractOperator,
                  eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

  lo2 = getBaseLO(lo)
  @assert typeof(lo2) <: PetscMatFreeLO
  lo2.ctx = (mesh, sbp, eqn, opts, lo, ctx_residual, t)
  MatShellSetContext(lo2.A, pointer_from_objref(lo2.ctx))

  return nothing
end

function applyLinearOperator_wrapper(A::PetscMat, x::PetscVec, b::PetscVec)

  ctx_ptr = MatShellGetContext(A)
  @assert ctx_ptr != C_NULL
  ctx = unsafe_pointer_to_objref(ctx_ptr)
  checkLOCtx(ctx)
  mesh = ctx[1]
  sbp = ctx[2]
  eqn = ctx[3]
  opts = ctx[4]
  lo = ctx[5]  # user defined LO
  ctx_residual = ctx[6]
  t = ctx[7]

  btmp = VecGetArray(b)
  xtmp = VecGetArrayRead(x)

  applyLinearOperator(lo, mesh, sbp, eqn, opts, ctx_residual, t, xtmp, btmp)

  VecRestoreArray(b, btmp)
  VecRestoreArrayRead(x, xtmp)

  return PetscErrorCode(0)
end

function applyLinearOperatorTranspose_wrapper(A::PetscMat, x::PetscVec,
                                              b::PetscVec)

  ctx_ptr = MatShellGetContext(A)
  @assert ctx_ptr != C_NULL
  ctx = unsafe_pointer_to_objref(ctx_ptr)
  checkLOCtx(ctx)
  mesh = ctx[1]
  sbp = ctx[2]
  eqn = ctx[3]
  opts = ctx[4]
  lo = ctx[5]  # user defined LO
  ctx_residual = ctx[6]
  t = ctx[7]

  btmp = VecGetArray(b)
  xtmp = VecGetArrayRead(x)

  applyLinearOperatorTranspose(lo, mesh, sbp, eqn, opts, ctx_residual, t,
                               xtmp, btmp)

  VecRestoreArray(b, btmp)
  VecRestoreArrayRead(x, xtmp)

  return PetscErrorCode(0)
end

"""
  This function verifies the [`PetscMatFreeLO`](@ref) ctx is valid, throwing
  an exception otherwise.

  Currently this only checks whether the members of the tuple have the right
  abstract supertypes

  **Inputs**

   * ctx: the ctx retrieved from the Petsc matrix shell object
"""
function checkLOCtx(ctx)

  mesh = ctx[1]
  sbp = ctx[2]
  eqn = ctx[3]
  opts = ctx[4]
  lo = ctx[5]  # user defined LO
  ctx_residual = ctx[6]
  t = ctx[7]

  lo2 = getBaseLO(lo)
  @assert typeof(mesh) <: AbstractMesh
  @assert typeof(sbp) <: AbstractOperator
  @assert typeof(eqn) <: AbstractSolutionData
  @assert typeof(lo) <: AbstractPetscMatFreeLO
  @assert typeof(lo2) <: PetscMatFreeLO
  @assert typeof(t) <: Number

  return nothing
end

function getBaseLO(lo::PetscMatFreeLO)

  # this is the bottom of the recursion tree
  return lo
end


text:
function transform_docstring(doc::AbstractString, label)
    output = IOBuffer()
    input = IOBuffer(doc)
    while !eof(input)
        ln = readline(input)
        if startswith(ln, "```julia")
            print(output, "```jldoctest ", label, "\n")
            isrepl = false
            while !eof(input)
                ln = readline(input)
                if startswith(ln, "```")
                    if !isrepl
                        print(
                            output,
                            """
                            nothing
                            # output
                            """,
                        )
                    end
                    print(output, ln, "\n")
                    break
                end
                print(output, ln, "\n")
                if startswith(ln, "julia> ")
                    isrepl = true
                end
            end
        else
            print(output, ln, "\n")
        end
    end
    return String(take!(output))
end

function define_docstrings()
    docstrings = [:Reagents => joinpath(dirname(@__DIR__), "README.md")]
    docsdir = joinpath(@__DIR__, "docs")
    for filename in readdir(docsdir)
        stem, ext = splitext(filename)
        ext == ".md" || continue
        name = Symbol(stem)
        name in names(Reagents, all = true) || continue
        push!(docstrings, name => joinpath(docsdir, filename))
    end
    n_auto_labels = 0
    for (name, path) in docstrings
        label = string(name)
        if match(r"^[a-z0-9_]+$"i, label) === nothing
            label = "reagent$n_auto_labels"
            n_auto_labels += 1
        end

        include_dependency(path)
        doc = read(path, String)
        doc = transform_docstring(doc, label)
        ex = :($Base.@doc $doc $name)
        ex.args[2]::LineNumberNode
        ex.args[2] = LineNumberNode(1, Symbol(path))
        Base.eval(Reagents, ex)
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetVMInstanceIDs <: SwaggerModel
    instanceIds::Any # spec type: Union{ Nothing, Vector{String} } # spec name: instanceIds

    function VirtualMachineScaleSetVMInstanceIDs(;instanceIds=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetVMInstanceIDs, Symbol("instanceIds"), instanceIds)
        setfield!(o, Symbol("instanceIds"), instanceIds)
        o
    end
end # type VirtualMachineScaleSetVMInstanceIDs

const _property_map_VirtualMachineScaleSetVMInstanceIDs = Dict{Symbol,Symbol}(Symbol("instanceIds")=>Symbol("instanceIds"))
const _property_types_VirtualMachineScaleSetVMInstanceIDs = Dict{Symbol,String}(Symbol("instanceIds")=>"Vector{String}")
Base.propertynames(::Type{ VirtualMachineScaleSetVMInstanceIDs }) = collect(keys(_property_map_VirtualMachineScaleSetVMInstanceIDs))
Swagger.property_type(::Type{ VirtualMachineScaleSetVMInstanceIDs }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetVMInstanceIDs[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetVMInstanceIDs }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetVMInstanceIDs[property_name]

function check_required(o::VirtualMachineScaleSetVMInstanceIDs)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetVMInstanceIDs }, name::Symbol, val)
end


text:
function dlpRead(fName::AbstractString)::tsParm
    o=tsParm()

    # A. Read dlp file - all into memory (unless really big file this
    #     should not be a problem when optimisation is done

    # A1. Get data to preallocate memory
    #     Also count the number of data lines
    fID::IOStream=open(fName)
    nLines::Int64=0
    ndLines::Int64=0
    iDataSt::Int64=0
    iDataEd::Int64=0
    foundData::Bool=false
    rLines::Array{AbstractString,1}=[]
    while !eof(fID)
      in::AbstractString=readline(fID)
      nLines+=1
      push!(rLines,in)
      if foundData
        m=match(r"^\%",in)    # Check this is not a comment line
        if m === nothing
          m=match(r"[0-9]+",in) # Check there is something on the line - not just blanks
          if m !== nothing
	    (ndLines==0) ? iDataSt=nLines : 0
            ndLines+=1
	  else
	    if iDataEd==0
	      (ndLines>0) ? iDataEd=nLines-1 : 0
	    end
	  end
        end
      end
      m=match(r"\[data\]"i,in)  # Check if this is the start of the data section
      if m !== nothing
        foundData=true
      end
    end
    (iDataEd==0) ? iDataEd=nLines : 0
    close(fID)
    println("Checked file $(fName):")
    println("  Found $(nLines) lines")
    println("  Found data lines from $(iDataSt) to $(iDataEd)")

    # B. Now make some interpretations of the Parameter sections
    function getParNm(inStr) m=match(r".*:(.*)]",inStr); a=m.captures[1]; end
    iP=find(x->contains(x,"[parameter"),rLines)
    parNms=Array(AbstractString,length(iP))
    for k=1:length(iP)
      parNms[k]=getParNm(rLines[iP[k]])
    end

    # C. Loop over each parameter section and gather metadata in a Dict
    function readParSect(iP::Int64, rLines::Array{AbstractString})
      c::Dict=Dict()
      lk=0
      while(!(chomp(rLines[iP+lk+1])==""))
        lk=lk+1
        m=match(r"(.* ) +(.*)",rLines[iP+lk])
        if m!==nothing
           a=replace(m.captures[1]," ","");
           b=replace(m.captures[2],r"\r","");
	   c[a]=b
        end
      end
      return c
    end

    local k::Int64=0
    o.Meta=Dict()
    for k=1:length(iP)
      o.Meta[parNms[k]]=readParSect(iP[k],rLines)
    end

    # C1. Set list of variables
    o.Name=Array(AbstractString,length(iP))
    i=0
    for k=1:length(iP)
      if haskey(o.Meta[parNms[k]],"column")
        i+=1
        o.Name[i]=getParNm(rLines[iP[k]])
      end
    end

    # D. Interpret in the dataset
    i=0
    df = Dates.DateFormat("dd-uuu-yyyy HH:MM");
    for k=iDataSt:iDataEd
       a=replace(rLines[k],"\r\n","")
       a=replace(a,","," ")
       kk=search(a,r" |,")
       kk=search(a,r" |,",kk[end]+1)
       tDate=a[1:kk[1]]
       a=replace(a[kk[1]:end],r" +"," ")
       a=replace(a,r"^ ","")
       a=replace(a,"\n","")
       jdd=readdlm(IOBuffer(a))
       if i==0
         o.sDate=[]
         o.Data=zeros(iDataEd-iDataSt+1,size(jdd,2))
       end
       i+=1
       o.Data[i,:]=jdd[1,:]
       push!(o.sDate,Dates.DateTime(tDate,df))
    end
    println("  Read $i data lines")
    return o

end

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
export Interpolator, interpolatePoint, interpolateMesh, gradientPoint, gradientPoint, gradientMesh
using Parameters
using StaticArrays
using ForwardDiff
include("./utils/RBFutils.jl")
#######################################################
#               Interpolate environment
#######################################################
"""
Interpolator environment setup :: mutable struct Interpolator
__________________________________________________________________________________________________
Parameters:                                                                                |
    X::Matrix           Matrix to set interpolation                                               |
    y::Vector           Vector with interpolate ready variables.                                  |
    ϵ::Float64          Epsilon to set basis function                    => Default : 1.          |
    basis_name::Strin   Selected for user ready basis                                             |
    φ::Function         Function of (distance between points ,  epsilon) => Default : phi_MQ(r,ϵ) |
    w::SVector          Static Vector containing the weights                                      |
__________________________________________________________________________________________________|
"""
@with_kw mutable struct Interpolator
    X::Matrix # Matrix to set interpolation
    y::Vector # Vector with interpolate ready variables
    ϵ::Float64 = 1.0 # Epsilon to set basis function
    basis_name::String = "MQ" #Selected for user ready basis
    φ::Function = _get_kernel(basis_name) # Function of (distance between points ,  epsilon)
    w::Vector = _build(X, y; φ = φ, ϵ = ϵ) # To setup the interpolation environment
end
#######################################################
#               Interpolate point
#######################################################
"""
Function to interpolate one point
__________________________________________________________________________________________________
Parameters:                                                                                |
    inter::Interpolator         An interpolation environment already setted up                    |
    p::Vector                   Coordinates of the point to interpolate                           |
Output:                                                                                |
    Interpolation results setted by y on Interpolator                                             |
__________________________________________________________________________________________________|
"""
function interpolatePoint(inter::Interpolator, p::Vector)
    ((dot(inter.w, [inter.φ(r(p, inter.X[i]), inter.ϵ) for i = 1:size(inter.X)[1]]) * std(inter.y))) + mean(inter.y)
end

#######################################################
#               Interpolate mesh
#######################################################

"""
Function to interpolate one point
__________________________________________________________________________________________________
Parameters:                                                                                |
    inter::Interpolator         An interpolation environment already setted up                    |
    ps::Matrix                  Matrix with the mesh to interpolate                               |
Output:                                                                                |
    Interpolation results setted by y on Interpolator                                             |
__________________________________________________________________________________________________|
"""

function interpolateMesh(inter::Interpolator, ps::Matrix)
    [interpolatePoint(inter, p[:]) for p in eachrow(ps)]
end
#######################################################
#               Interpolate gradient
#######################################################
"""
Function to interpolate one point's gradient
__________________________________________________________________________________________________
Parameters:                                                                                |
    inter::Interpolator         An interpolation environment already setted up                    |
    p::Vector                   Coordinates of the point to interpolate                           |
Output:                                                                                |
    Interpolation results setted by y on Interpolator and using ForwardDiff                       |
__________________________________________________________________________________________________|
"""
function gradientPoint(inter, p; N = 10000)
    f(r) = interpolatePoint(inter, r)
    ForwardDiff.gradient(f, p)
end


text: include("dp.jl")