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

module QuantumTomography

import Distributions.fit

export fit,
       #fitA,
       #fitB,
       predict,
       FreeLSStateTomo,
       LSStateTomo,
       MLStateTomo

using Convex, Distributions, SCS, QuantumInfo

function build_state_predictor(obs::Vector{Matrix})
    return reduce(vcat,[vec(o)' for o in obs])
end

function build_process_predictor(obs::Vector{Matrix}, prep::Vector{Matrix})
    exps = Matrix[ choi_liou_involution(vec(o)*vec(p)') for o in obs, p in prep ]
    return reduce(vcat, map(m->vec(m)', vec(exps)) )
end

"""
`FreeLSStateTomo`

Free (unconstrained) least-squares state tomography algorithm. It is
constructed from a vector of observables corresponding to
measurements that are performed on the state being reconstructed.
"""
type FreeLSStateTomo
    inputdim::Int
    outputdim::Int
    pred::Matrix
    function FreeLSStateTomo(obs::Vector)
        pred = build_state_predictor(obs)
        outputdim = size(pred,1)
        inputdim = size(pred,2)
        new(inputdim,outputdim,pred)
    end
end

"""
`predict`

Predict outcomes of a tomography experiment for a given state (density matrix).
"""
function predict(method::FreeLSStateTomo, state)
    method.pred*vec(state)
end

"""
`fit(method, ...)`

Reconstruct a state from observations (i.e., perform state tomography). Exactly which fitting
routine is used depends on the type of the `methos`. The possible types are

   + FreeLSStateTomo : unconstrained least-squares state tomography
   + LSStateTomo : least-squares state tomography constrained to physical states
   + MLStateTomo : maximum-likelihood, maximum-entropy state tomography constrained to physical states

"""
function fit(method::FreeLSStateTomo, means::Vector{Float64}, vars::Vector{Float64}=-ones(length(means)); algorithm=:OLS)
    if length(means) != method.outputdim
        error("The number of expected means does not match the required number of experiments")
    end
    d = round(Int, method.inputdim |> sqrt)
    if algorithm==:OLS
        reg = method.pred\means
        return reshape(reg,d,d), norm(method.pred*reg-means,2)/length(means), :Optimal
    elseif algorithm==:GLS
        if any(vars<0)
            error("Variances must be positive for generalized least squares.")
        end
        reg = (sqrt(vars)\method.pred)\means
        return reshape((sqrt(vars)\method.pred)\means,d,d), sqrt(dot(method.pred*reg-means,diagm(vars)\(method.pred*reg-means)))/length(means), :Optimal
    else
        error("Unrecognized method for least squares state tomography")
    end
end

"""
`LSStateTomo`

Type corresponding to constrained least-squares state tomography
algorithm. It is constructed from a vector of observables
corresponding to measurements that are performed on the state being
reconstructed.

The outcome of said measurements for a given density matrix ρ can be
computed using the `predict` function , while the density matrix can
be estimated from observations (including variances) by using the
`fit` function.

"""
type LSStateTomo
    inputdim::Int
    outputdim::Int
    realpred::Matrix
    function LSStateTomo(obs::Vector)
        pred = build_state_predictor(obs)
        outputdim = size(pred,1)
        inputdim = size(pred,2)
        realpred = [real(pred) imag(pred)];
        new(inputdim,outputdim,realpred)
    end
end

function predict(method::LSStateTomo, state)
    (method.realpred[:,1:round(Int,end/2)]+1im*method.realpred[:,round(Int,end/2)+1:end])*vec(state)
end

function fit(method::LSStateTomo,
             means::Vector{Float64},
             vars::Vector{Float64};
             solver = SCSSolver(verbose=0, max_iters=10_000, eps = 1e-8))

    if length(means) != length(vars) || method.outputdim != length(means)
        error("Size of observations and/or predictons do not match.")
    end
    dsq = method.inputdim
    d = round(Int,sqrt(dsq))

    # We assume that the predictions are always real-valued
    # and we need to do the complex->real translation manually since
    # Convex.jl does not support complex numbers yet
    ivars = 1./sqrt(vars)

    ρr = Variable(d,d)
    ρi = Variable(d,d)

    constraints = trace(ρr) == 1
    constraints += trace(ρi) == 0
    constraints += isposdef([ρr ρi; -ρi ρr])

    # TODO: use quad_form instead of vecnorm? Have 1/vars are diagonal quadratic form
    problem = minimize( vecnorm( (means - method.realpred*[vec(ρr); vec(ρi)]) .* ivars, 2), constraints )

    solve!(problem, solver)

    return (ρr.value - 1im*ρi.value), problem.optval^2, problem.status
end

"""
`MLStateTomo`

Maximum-likelihood maximum-entropy quantum state tomography algorithm. It is
constructed from a collection of observables corresponding to
measurements that are performed on the state being reconstructed, as
well as a hedging factor β. If β=0, no hedging is applied. If β > 0
either a log determinant or log minimum eigenvalue penalty is applied.

**AT THE MOMENT HEDING IS ONLY APPLIED IN fitA, AND IS NOT CURRENTLY WORKING**
"""
type MLStateTomo
    effects::Vector
    dim::Int64
    β::Float64
    function MLStateTomo(v::Vector,β=0.0)
        for e in v
            if !ishermitian(e) || !ispossemidef(e) || real(trace(e))>1
                error("MLStateTomo state tomography is parameterized by POVM effects only.")
            end
        end
        if !all([size(e,1)==size(e,2) for e in v]) || !all([size(v[1],1)==size(e,1) for e in v])
                error("All effects must be square matrices, and they must have have the same dimension.")
        end
        if β < 0
            error("Hedging penalty must be positive.")
        end
        new(v,size(v[1],1),β)
    end
end

function predict(method::MLStateTomo,ρ)
    Float64[ real(trace(ρ*e)) for e in method.effects]
end

"""
`fitA(method::MLStateTomo, freq, solver=SCSSolver)`

This is a state tomography fitting routine using convex optimization. It's use is currently
discouraged simply because it is much slower than the iterative solver, and often does not
conver to a solution that meets the optimality criteria.
"""
function fitA(method::MLStateTomo,
             freq::Vector;
             solver = SCSSolver(verbose=0, max_iters=100_000, eps = 1e-8))

    if length(method.effects) != length(freq)
        error("Vector of counts and vector of effects must have same length, but length(counts) == $(length(counts)) != $(length(method.effects))")
    end

    d = method.dim

    # ρr = Variable(d,d)
    ρr = Semidefinite(d)
    ρi = Variable(d,d)

    ϕ(m) = [real(m) -imag(m); imag(m) real(m)];
    ϕ(r,i) = [r i; -i r]

    obj = freq[1] * log(trace([ρr ρi; -ρi ρr]*ϕ(method.effects[1])))
    for i=2:length(method.effects)
        obj += freq[i] * log(trace([ρr ρi; -ρi ρr]*ϕ(method.effects[i])))
    end
    # add hedging regularization
    if method.β != 0.0
        obj += method.β*log(lambdamin([ρr ρi; -ρi ρr]))
    end

    constraints = trace(ρr) == 1
    constraints += trace(ρi) == 0
    constraints += isposdef([ρr ρi; -ρi ρr])

    problem = maximize(obj, constraints)

    solve!(problem, solver)

    return ρr.value-im*ρi.value, problem.optval, problem.status
end

"""
`fitB(method::MLStateTomo, freq, solver=SCSSolver)`

This is a state tomography fitting routine using convex optimization. It's use is currently
discouraged simply because it is much slower than the iterative solver, and often does not
conver to a solution that meets the optimality criteria, much like fitB.
"""
function fitB(method::MLStateTomo,
             freq::Vector;
             solver = SCSSolver(verbose=0, max_iters=100_000, eps = 1e-8))

    if length(method.effects) != length(freq)
        error("Vector of counts and vector of effects must have same length, but length(counts) == $(length(counts)) != $(length(method.effects))")
    end

    d = method.dim

    ρr = Semidefinite(d) #Variable(d,d)
    ρi = Variable(d,d)
    p  = Variable(length(method.effects),Positive())

    ϕ(m) = [real(m) -imag(m); imag(m) real(m)];
    ϕ(r,i) = [r i; -i r]

    obj = freq[1] * log(p[1])
    for i=2:length(method.effects)
        obj += freq[i] * log(p[i])
    end
    #obj += method.β!=0.0 ? method.β*logdet([ρr ρi; -ρi ρr]) : 0.0
    #obj += method.β!=0.0 ? method.β*log(lambdamin([ρr ρi; -ρi ρr])) : 0.0

    constraints = trace(ρr) == 1
    constraints += trace(ρi) == 0
    for i=1:length(method.effects)
        constraints += p[i] == trace([ρr ρi; -ρi ρr]*ϕ(method.effects[i]))/2
    end
    constraints += isposdef([ρr ρi; -ρi ρr])

    problem = maximize(obj, constraints)

    solve!(problem, solver)

    return ρr.value-im*ρi.value, problem.optval, problem.status
end

function R(ρ,E,obs)
    pr = Float64[real(trace(ρ*Ei)) for Ei in E]
    R = obs[1]/pr[1]*E[1]
    for i in 2:length(E)
        R += obs[i]/pr[i]*E[i]
    end
    R
end

function LL(ρ,E,obs)
    pr = Float64[real(trace(ρ*Ei)) for Ei in E]
    sum(obs.*log(pr))
end

function fit(method::MLStateTomo,
             freq::Vector;
             δ=.005, # a more natural parameterization of "dilution"
             λ=0.0,  # entropy penalty
             tol=1e-9,
             maxiter=100_000,
             ρ0 = eye(Complex128,method.dim))
    ϵ=1/δ
    iter = 1
    ρk = copy(ρ0)
    ρktemp = similar(ρk)
    ρkm = similar(ρ0) #Array(Complex128,method.dim,method.dim)
    status = :Optimal
    while true
        copy!(ρkm,ρk)
        if iter >= maxiter
            status = :MaxIter
            break
        end
        # diluted iterative scheme for ML (from PRA 75 042108 2007)
        Tk = λ > 0.0 ?
             (1+ϵ*(R(ρk,method.effects,freq)-eye(method.dim)-λ*(logm(ρk)-trace(ρk*logm(ρk)))))/(1+ϵ) :
             (1+ϵ*R(ρk,method.effects,freq))/(1+ϵ)
        A_mul_B!(ρktemp,ρk,(1+ϵ*Tk)/(1+ϵ)) # ρk = ρk * (1+ϵRk)/(1+ϵ)
        A_mul_B!(ρk,(1+ϵ*Tk)/(1+ϵ),ρktemp) # ρk = (1+ϵRk) * ρk/(1+ϵ)
        trnormalize!(ρk)
        if vecnorm(ρk-ρkm)/method.dim^2 < tol
            status = :Optimal
            break
        end
        iter += 1
    end
    return ρk, LL(ρk,method.effects,freq), status
end

function trb_sop(da,db)
    sop = spzeros(da^2,(da*db)^2)
    for i=0:da-1
        for j=0:da-1
            for k=0:db-1
                sop += vec(ketbra(i,j,da))*vec(kron(ketbra(i,j,da),ketbra(k,k,db)))'
            end
        end
    end
    return sop
end

# TODO: what about the (unobservable) traceful component
function qpt_lsq(pred::Matrix, means::Vector{Float64}, vars::Vector{Float64}; method=:OLS)
    d = Int(shape(pred,1) |> sqrt |> round)
    if method==:OLS
        return reshape(pred\means,d,d)
    elseif method==:GLS
        return reshape((sqrt(vars)\pred)\means,d,d)
    else
        error("Unrecognized method for least squares process tomography")
    end
end

type LSProcessTomo
    inputdim::Int
    outputdim::Int
    realpred::Vector{AbstractMatrix}
end

# For QPT, we write the predictor as operating on Choi-Jamilokoski
# matrices.  This is a bit awkward in comparisson to using the
# Liouville/natural representation, but it gets around some of the
# limitations of Convex.jl, and it is also much more efficient.
function qpt_ml(pred::Matrix, means::Vector{Float64}, vars::Vector{Float64})
    if length(means) != length(vars) || size(pred,1) != length(means)
        error("Size of observations and/or predictons do not match.")
    end
    d4 = size(pred,2)
    d2 = Int(sqrt(d4))
    d  = Int(sqrt(d2))

    # We assume that the predictions are always real-valued
    # and we need to do the complex->real translation manually since
    # Convex.jl does not support complex numbers yet
    rpred = [real(pred) imag(pred)];
    ivars = 1./sqrt(vars)

    ptrb = trb_sop(d,d)

    ρr = Variable(d2,d2)
    ρi = Variable(d2,d2)

    problem = minimize( vecnorm( (means - rpred*[vec(ρr); vec(ρi)]) .* ivars, 2 )^2 )

    problem.constraints += isposdef([ρr ρi; -ρi ρr])
    problem.constraints += trace(ρi) == 0
    problem.constraints += reshape(ptrb*vec(ρr),d,d) == eye(d)
    problem.constraints += reshape(ptrb*vec(ρi),d,d) == zeros(d,d)

    solve!(problem, SCSSolver(verbose=0))

    return (ρr.value - 1im*ρi.value), problem.optval, problem.status
end

function qst_ml(obs::Vector{Matrix}, means::Vector{Float64}, vars::Vector{Float64})
    #pred = reduce(vcat,[vec(o)' for o in obs])
    pred = build_state_predictor(obs)
    return qst_ml(pred, means, vars)
end

function qpt_ml(obs::Vector{Matrix}, states::Vector{Matrix}, means::Vector{Float64}, vars::Vector{Float64})
    #exps = Matrix[ choi_liou_involution(vec(o)*vec(p)') for o in obs, p in prep ]
    #pred = reduce(vcat, map(m->vec(m)', vec(exps)) )
    pred = build_process_predictor(obs, prep)
    return qpt_ml(pred, meas, vars)
end

end # module


text:
@doc raw"""
    LinearScalarAdvectionEquation1D

The linear scalar advection equation
```math
\partial_t u + a \partial_1 u  = 0
```
in one space dimension with constant velocity `a`.
"""
struct LinearScalarAdvectionEquation1D <: AbstractLinearScalarAdvectionEquation{1, 1}
  sources::String
  advectionvelocity::SVector{1, Float64}
end

function LinearScalarAdvectionEquation1D()
  sources = parameter("sources", "none")
  a = convert(SVector{1,Float64}, parameter("advectionvelocity"))
  LinearScalarAdvectionEquation1D(sources, a)
end

get_name(::LinearScalarAdvectionEquation1D) = "LinearScalarAdvectionEquation1D"
varnames_cons(::LinearScalarAdvectionEquation1D) = SVector("scalar")
varnames_prim(::LinearScalarAdvectionEquation1D) = SVector("scalar")

# Set initial conditions at physical location `x` for time `t`
function initial_conditions_gauss(x, t, equation::LinearScalarAdvectionEquation1D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  return @SVector [exp(-(x_trans[1]^2 ))]
end

function initial_conditions_convergence_test(x, t, equation::LinearScalarAdvectionEquation1D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  c = 1.0
  A = 0.5
  L = 2
  f = 1/L
  omega = 2 * pi * f
  scalar = c + A * sin(omega * sum(x_trans))
  return @SVector [scalar]
end

function initial_conditions_sin(x, t, equation::LinearScalarAdvectionEquation1D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  scalar = sin(2 * pi * x_trans[1])
  return @SVector [scalar]
end

function initial_conditions_constant(x, t, equation::LinearScalarAdvectionEquation1D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  return @SVector [2.0]
end

function initial_conditions_linear_x(x, t, equation::LinearScalarAdvectionEquation1D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  return @SVector [x_trans[1]]
end

# Apply boundary conditions
function boundary_conditions_linear_x(u_inner, orientation, direction, x, t, surface_flux_function,
                                      equation::LinearScalarAdvectionEquation1D)
  u_boundary = initial_conditions_linear_x(x, t, equation)

  # Calculate boundary flux
  if direction == 2  # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

function boundary_conditions_gauss(u_inner, orientation, direction, x, t, surface_flux_function,
                                   equation::LinearScalarAdvectionEquation1D)
  u_boundary = initial_conditions_gauss(x, t, equation)

  # Calculate boundary flux
  if direction == 2  # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

function boundary_conditions_convergence_test(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equation::LinearScalarAdvectionEquation1D)
  u_boundary = initial_conditions_convergence_test(x, t, equation)

  # Calculate boundary flux
  if direction == 2  # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

# Pre-defined source terms should be implemented as
# function source_terms_WHATEVER(ut, u, x, element_id, t, n_nodes, equation::LinearScalarAdvectionEquation2D)

# Calculate 1D flux in for a single point
@inline function calcflux(u, orientation, equation::LinearScalarAdvectionEquation1D)
  a = equation.advectionvelocity[orientation]
  return a * u
end

function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::LinearScalarAdvectionEquation1D)
  a = equation.advectionvelocity[orientation]
  return 0.5 * ( a * (u_ll + u_rr) - abs(a) * (u_rr - u_ll) )
end

# Determine maximum stable time step based on polynomial degree and CFL number
function calc_max_dt(u, element_id, invjacobian, cfl,
                     equation::LinearScalarAdvectionEquation1D, dg)
  λ_max = maximum(abs, equation.advectionvelocity)
  return cfl * 2 / (nnodes(dg) * invjacobian * λ_max)
end

# Convert conservative variables to primitive
@inline cons2prim(u, equation::LinearScalarAdvectionEquation1D) = u

# Convert conservative variables to entropy variables
@inline cons2entropy(u, equation::LinearScalarAdvectionEquation1D) = u

# Calculate entropy for a conservative state `cons`
@inline entropy(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5 * u^2
@inline entropy(u, equation::LinearScalarAdvectionEquation1D) = entropy(u[1], equation)

# Calculate total energy for a conservative state `cons`
@inline energy_total(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5 * u^2
@inline energy_total(u, equation::LinearScalarAdvectionEquation1D) = energy_total(u[1], equation)

ext: jl


text:
module AMLPipelineR

using AMLPipelineBase
using AMLPipelineBase: fit!,transform!,fit_transform!
using AMLPipelineBase: AbsTypes, Utils, BaselineModels, Pipelines
using AMLPipelineBase: Machine, Learner, Transformer, Workflow, Computer
using AMLPipelineBase: BaseFilters, FeatureSelectors, DecisionTreeLearners
using AMLPipelineBase: EnsembleMethods, CrossValidators
using AMLPipelineBase: NARemovers

import AMLPipelineBase.AbsTypes: fit!, transform!

export fit!, transform!, fit_transform!
export Machine, Learner, Transformer, Workflow, Computer
export Baseline, Identity
export Imputer,OneHotEncoder,Wrapper
export PrunedTree,RandomForest,Adaboost
export VoteEnsemble, StackEnsemble, BestLearner
export FeatureSelector, CatFeatureSelector, NumFeatureSelector
export CatNumDiscriminator
export crossvalidate
export NARemover
export @pipeline @pipelinex, @pipelinez
export Pipeline, ComboPipeline

include("rmachine.jl")
using .RMachine
export CRTLearner

end #module

max_stars_repo_name: ppalmes/AMLPipelineR

max_issues_repo_name: ppalmes/AMLPipelineR

max_forks_repo_name: ppalmes/AMLPipelineR


text:
mutable struct PcapStreamWriter{Dst <: IO} <: PcapWriter
    dst::Dst

    function PcapStreamWriter{Dst}(dst::Dst; thiszone = 0, snaplen = 65535) where {Dst <: IO}
        h = PcapHeader(
            0xa1b23c4d,
            0x0002,
            0x0004,
            thiszone,
            0,
            snaplen,
            LINKTYPE_ETHERNET)
        write(dst, Ref(h))
        new(dst)
    end
end

PcapStreamWriter(io::IO; kwargs...) = PcapStreamWriter{typeof(io)}(io; kwargs...)
PcapStreamWriter(path::AbstractString; kwargs...) = PcapStreamWriter(open(path, "w"); kwargs...)

Base.close(x::PcapStreamWriter) = close(x.dst)

function Base.write(x::PcapStreamWriter, timestamp::UnixTime, data::DenseArray{UInt8})
    sec, nsec = fldmod(Dates.value(timestamp), 1_000_000_000)
    data_length = length(data)
    h = RecordHeader(sec, nsec, data_length, data_length)
    write(x.dst, Ref(h))
    write(x.dst, data)
    nothing
end


text:
using Cxx

@static if is_apple()
    if @__DIR__ != nothing
        const DLIB_PATH = "$(@__DIR__)/../deps/build/dlib/dlib/build/libdlib.dylib"
    else
        const DLIB_PATH = Pkg.dir("Dlib", "deps", "build", "dlib", "dlib", "build", "libdlib.dylib")
    end
end

function init()
    Libdl.dlopen(DLIB_PATH, Libdl.RTLD_GLOBAL)
end


text:
# Deterministic equivalent generation #
# ========================== #
function generate!(stochasticprogram::StochasticProgram{N}, structure::DeterministicEquivalent{N}) where N
    # Set the optimizer
    structure.model.moi_backend = optimizer(stochasticprogram)
    # Prepare decisions
    structure.model.ext[:stage_map] = Dict{MOI.VariableIndex, Int}()
    structure.model.ext[:decisions] = structure.decisions
    add_decision_bridges!(structure.model)
    # Generate all stages
    for stage in 1:N
        generate!(stochasticprogram, structure, stage)
    end
    return nothing
end

function generate!(stochasticprogram::StochasticProgram{N}, structure::DeterministicEquivalent{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    dep_model = structure.model
    if stage == 1
        # Define first-stage problem
        generator(stochasticprogram, :stage_1)(dep_model, stage_parameters(stochasticprogram, 1))
        # Cache first-stage objective
        obj = objective_function(dep_model)
        sense = objective_sense(dep_model)
        push!(structure.sub_objectives[1], (sense, moi_function(obj)))
        # Bookkeeping
        for (objkey, obj) in object_dictionary(dep_model)
            if isa(obj, DecisionRef)
                # Get common index in proxy and add to variable map
                proxy_var = proxy(stochasticprogram, stage)[objkey]
                for s in 2:N
                    # TODO: Update for multi-stage
                    for i in 1:num_scenarios(structure, s)
                        structure.variable_map[(index(proxy_var), i)] = index(obj)
                    end
                end
            elseif isa(obj, AbstractArray{<:DecisionRef})
                # Get common index in proxy
                proxy_var = proxy(stochasticprogram, stage)[objkey]
                for (var, proxy) in zip(obj, proxy_var)
                    # Update variable map
                    for s in 2:N
                        # TODO: Update for multi-stage
                        for i in 1:num_scenarios(structure, s)
                            structure.variable_map[(index(proxy), i)] = index(var)
                        end
                    end
                end
            end
        end
    else
        # Sanity check on scenario probabilities
        if num_scenarios(stochasticprogram, stage) > 0
            p = stage_probability(stochasticprogram, stage)
            abs(p - 1.0) <= 1e-6 || @warn "Scenario probabilities do not add up to one. The probability sum is given by $p"
        end
        # Cache current objective and sense
        dep_obj = objective_function(dep_model)
        obj_sense = objective_sense(dep_model)
        obj_sense = obj_sense == MOI.FEASIBILITY_SENSE ? MOI.MIN_SENSE : obj_sense
        # Null objective temporarily in case subproblem objectives are zero
        @objective(dep_model, obj_sense, 0)
        # Define second-stage problems, renaming variables according to scenario.
        stage_two_params = stage_parameters(stochasticprogram, 2)
        visited_objs = collect(keys(object_dictionary(dep_model)))
        seen_constraints = CI[]
        # Do not need to map any first-stage decision constraints
        for (F,S) in MOI.get(structure.model, MOI.ListOfConstraints())
            if is_decision_type(F)
                append!(seen_constraints, MOI.get(structure.model, MOI.ListOfConstraintIndices{F,S}()))
            end
        end
        # Loop through scenarios and incrementally build deterministic equivalent
        for (i, scenario) in enumerate(scenarios(stochasticprogram))
            # Generate model information for scenario i
            stage_key = Symbol(:stage_, stage)
            generator(stochasticprogram, stage_key)(dep_model, stage_parameters(stochasticprogram, stage), scenario)
            # Update objective and cache the subobjective function for scenario i
            sub_sense = objective_sense(dep_model)
            sub_obj = objective_function(dep_model)
            if obj_sense == sub_sense
                dep_obj += probability(scenario) * objective_function(dep_model)
            else
                dep_obj -= probability(scenario) * objective_function(dep_model)
            end
            push!(structure.sub_objectives[stage], (sub_sense, moi_function(sub_obj)))
            # Bookkeeping for objects added in scenario i
            for (objkey,obj) in filter(kv->kv.first ∉ visited_objs, object_dictionary(dep_model))
                newkey = if isa(obj, VariableRef) || isa(obj, DecisionRef)
                    if isa(obj, DecisionRef)
                        # Get common index in proxy and add to variable map
                        proxy_var = proxy(stochasticprogram, stage)[objkey]
                        structure.variable_map[(index(proxy_var), i)] = index(obj)
                    end
                    # Update variable name to reflect stage and scenario
                    varname = if N > 2
                        varname = add_subscript(add_subscript(JuMP.name(obj), stage), i)
                    else
                        varname = add_subscript(JuMP.name(obj), i)
                    end
                    set_name(obj, varname)
                    # Return new key
                    newkey = Symbol(varname)
                elseif isa(obj, AbstractArray{<:VariableRef}) || isa(obj, AbstractArray{<:DecisionRef})
                    # Update object name to reflect stage and scenario
                    arrayname = if N > 2
                        arrayname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        arrayname = add_subscript(objkey, i)
                    end
                    # Get common indices from proxy
                    proxy_var = proxy(stochasticprogram, stage)[objkey]
                    # Handle each individual variable in array
                    for (var, proxy) in zip(obj, proxy_var)
                        # Update variable map
                        if isa(obj, AbstractArray{<:DecisionRef})
                            structure.variable_map[(index(proxy), i)] = index(var)
                        end
                        # Update variable name to reflect stage and scenario
                        splitname = split(name(var), "[")
                        varname = if N > 2
                            varname = add_subscript(splitname[1], stage)
                        else
                            varname = splitname[1]
                        end
                        set_name(var, @sprintf("%s[%s", add_subscript(varname,i), splitname[2]))
                    end
                    # Return new key
                    newkey = Symbol(arrayname)
                elseif isa(obj, JuMP.ConstraintRef) ||
                       isa(obj, AbstractArray{<:ConstraintRef}) ||
                       isa(obj, AbstractArray{<:GenericAffExpr}) ||
                       isa(obj, AbstractArray{<:DecisionAffExpr})
                    # Update obj name to reflect stage and scenario
                    arrayname = if N > 2
                        arrayname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        arrayname = add_subscript(objkey, i)
                    end
                    # Return new key
                    newkey = Symbol(arrayname)
                elseif isa(obj, AbstractJuMPScalar)
                    # Update obj name to reflect stage and scenario
                    exprname = if N > 2
                        exprname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        exprname = add_subscript(objkey, i)
                    end
                    # Return new key
                    newkey = Symbol(exprname)
                else
                    continue
                end
                # Change name of obj to newkey to avoid collisions
                dep_model.obj_dict[newkey] = obj
                delete!(dep_model.obj_dict, objkey)
                # Bookkeep newkey to avoid handling again
                push!(visited_objs, newkey)
            end
            # Update constraint map
            for (F,S) in MOI.get(proxy(stochasticprogram, stage), MOI.ListOfConstraints())
                if is_decision_type(F)
                    constraints =  filter(MOI.get(structure.model, MOI.ListOfConstraintIndices{F,S}())) do ci
                        !(ci in seen_constraints)
                    end
                    proxy_constraints = MOI.get(proxy(stochasticprogram, stage), MOI.ListOfConstraintIndices{F,S}())
                    for (proxy,ci) in zip(proxy_constraints, constraints)
                        structure.constraint_map[(proxy, i)] = typeof(ci)(ci.value)
                        push!(seen_constraints, ci)
                    end
                end
            end
        end
        # Update objective
        set_objective_function(dep_model, dep_obj)
        set_objective_sense(dep_model, obj_sense)
    end
    return nothing
end

function clear!(dep::DeterministicEquivalent)
    # Clear decisions
    map(clear!, dep.decisions)
    # Clear subobjectives
    map(empty!, dep.sub_objectives)
    # Clear deterministic equivalent model
    empty!(dep.model)
    return nothing
end


text:
using RoboDojo
using LinearAlgebra
using DirectTrajectoryOptimization
const DTO = DirectTrajectoryOptimization

function hopper_dyn(mass_matrix, dynamics_bias, h, y, x, u, w)
    model = RoboDojo.hopper

    # dimensions
    nq = model.nq
    nu = model.nu

    # configurations

    q1⁻ = x[1:nq]
    q2⁻ = x[nq .+ (1:nq)]
    q2⁺ = y[1:nq]
    q3⁺ = y[nq .+ (1:nq)]

    # control
    u_control = u[1:nu]
    γ = u[nu .+ (1:4)]
    β = u[nu + 4 .+ (1:4)]

    E = [1.0 -1.0] # friction mapping
    J = RoboDojo.contact_jacobian(model, q2⁺)
    λ = transpose(J) * [[E * β[1:2]; γ[1]];
                        [E * β[3:4]; γ[2]];
                         γ[3:4]]
    λ[3] += (model.body_radius * E * β[1:2])[1] # friction on body creates a moment

    [
     q2⁺ - q2⁻;
     RoboDojo.dynamics(model, mass_matrix, dynamics_bias,
        h, q1⁻, q2⁺, u_control, zeros(model.nw), λ, q3⁺)
    ]
end

function hopper_dyn1(mass_matrix, dynamics_bias, h, y, x, u, w)
    nx = 8
    [
     hopper_dyn(mass_matrix, dynamics_bias, h, y, x, u, w);
     y[nx .+ (1:5)] - [u[2 .+ (1:4)]; u[end]];
     y[nx + 5 .+ (1:nx)] - x
    ]
end

function hopper_dynt(mass_matrix, dynamics_bias, h, y, x, u, w)
    nx = 8
    [
     hopper_dyn(mass_matrix, dynamics_bias, h, y, x, u, w);
     y[nx .+ (1:5)] - [u[2 .+ (1:4)]; u[end]];
     y[nx + 5 .+ (1:nx)] - x[nx + 5 .+ (1:nx)]
    ]
end

function contact_constraints_inequality_1(h, x, u, w)
    model = RoboDojo.hopper

    nq = model.nq
    nu = model.nu
    nx = 2nq

    q2 = x[1:nq]
    q3 = x[nq .+ (1:nq)]

    u_control = u[1:nu]
    γ = u[nu .+ (1:4)]
    β = u[nu + 4 .+ (1:4)]
    ψ = u[nu + 4 + 4 .+ (1:2)]
    η = u[nu + 4 + 4 + 2 .+ (1:4)]
    sα = u[nu + 4 + 4 + 2 + 4 .+ (1:1)]

    ϕ = RoboDojo.signed_distance(model, q3)

    μ = [model.friction_body_world; model.friction_foot_world]
    fc = μ .* γ[1:2] - [sum(β[1:2]); sum(β[3:4])]

    [
     -ϕ;
     -fc;
     β .* η .- sα;
     ψ .* fc  .- sα;
    ]
end

function contact_constraints_inequality_t(h, x, u, w)
    model = RoboDojo.hopper

    nq = model.nq
    nu = model.nu
    nx = 2nq

    q2 = x[1:nq]
    q3 = x[nq .+ (1:nq)]

    γ = u[nu .+ (1:4)]
    β = u[nu + 4 .+ (1:4)]
    ψ = u[nu + 4 + 4 .+ (1:2)]
    η = u[nu + 4 + 4 + 2 .+ (1:4)]
    sα = u[nu + 4 + 4 + 2 + 4 .+ (1:1)]

    ϕ = RoboDojo.signed_distance(model, q3)
    γ⁻ = x[nx .+ (1:4)]
    sα⁻ = x[nx + 4 .+ (1:1)]

    μ = [model.friction_body_world; model.friction_foot_world]
    fc = μ .* γ[1:2] - [sum(β[1:2]); sum(β[3:4])]

    [
     -ϕ;
     -fc;
     γ⁻ .* ϕ .- sα⁻;
     β .* η .- sα;
     ψ .* fc  .- sα;
    ]
end

function contact_constraints_inequality_T(h, x, u, w)
    model = RoboDojo.hopper

    nq = model.nq
    nx = 2nq

    q2 = x[1:nq]
    q3 = x[nq .+ (1:nq)]

    ϕ = RoboDojo.signed_distance(model, q3)
    γ⁻ = x[nx .+ (1:4)]
    sα⁻ = x[nx + 4 .+ (1:1)]

    [
     -ϕ;
     γ⁻ .* ϕ .- sα⁻;
    ]
end

function contact_constraints_equality(h, x, u, w)
    model = RoboDojo.hopper

    nq = model.nq
    nu = model.nu

    q2 = x[1:nq]
    q3 = x[nq .+ (1:nq)]

    γ = u[nu .+ (1:4)]
    β = u[nu + 4 .+ (1:4)]
    ψ = u[nu + 4 + 4 .+ (1:2)]
    η = u[nu + 4 + 4 + 2 .+ (1:4)]

    v = (q3 - q2) ./ h[1]
    vT_body = v[1] + model.body_radius * v[3]
    vT_foot = (RoboDojo.kinematics_foot_jacobian(model, q3) * v)[1]
    vT = [vT_body; -vT_body; vT_foot; -vT_foot]

    ψ_stack = [ψ[1] * ones(2); ψ[2] * ones(2)]

    [
     η - vT - ψ_stack;
    ]
end

# ## horizon
T = 21
h = 0.05

# ## hopper
nx = 2 * RoboDojo.hopper.nq
nu = RoboDojo.hopper.nu + 4 + 4 + 2 + 4 + 1

# ## model
mass_matrix, dynamics_bias = RoboDojo.codegen_dynamics(RoboDojo.hopper)
d1 = DTO.Dynamics((y, x, u, w) -> hopper_dyn1(mass_matrix, dynamics_bias, [h], y, x, u, w), 2 * nx + 5, nx, nu)
dt = DTO.Dynamics((y, x, u, w) -> hopper_dynt(mass_matrix, dynamics_bias, [h], y, x, u, w), 2 * nx + 5, 2 * nx + 5, nu)

dyn = [d1, [dt for t = 2:T-1]...];

# ## initial conditions
q1 = [0.0; 0.5 + RoboDojo.hopper.foot_radius; 0.0; 0.5]
qM = [0.5; 0.5 + RoboDojo.hopper.foot_radius; 0.0; 0.5]
qT = [1.0; 0.5 + RoboDojo.hopper.foot_radius; 0.0; 0.5]
q_ref = [0.5; 0.75 + RoboDojo.hopper.foot_radius; 0.0; 0.25]

x1 = [q1; q1]
xM = [qM; qM]
xT = [qT; qT]
x_ref = [q_ref; q_ref]

# ## gate
GAIT = 1
GAIT = 2
GAIT = 3

if GAIT == 1
	r_cost = 1.0e-1
	q_cost = 1.0e-1
elseif GAIT == 2
	r_cost = 1.0
	q_cost = 1.0
elseif GAIT == 3
	r_cost = 1.0e-3
	q_cost = 1.0e-1
end

# ## objective
function obj1(x, u, w)
	J = 0.0
	J += 0.5 * transpose(x - x_ref) * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0]) * (x - x_ref)
	J += 0.5 * transpose(u) * Diagonal([r_cost * ones(RoboDojo.hopper.nu); zeros(nu - RoboDojo.hopper.nu)]) * u
    J += 1000.0 * u[nu]
	return J
end

function objt(x, u, w)
	J = 0.0
	J += 0.5 * transpose(x[1:nx] - x_ref) * Diagonal(q_cost * [1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0]) * (x[1:nx] - x_ref)
	J += 0.5 * transpose(u) * Diagonal([r_cost * ones(RoboDojo.hopper.nu); zeros(nu - RoboDojo.hopper.nu)]) * u
    J += 1000.0 * u[nu]
	return J
end

function objT(x, u, w)
	J = 0.0
	J += 0.5 * transpose(x[1:nx] - x_ref) * Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]) * (x[1:nx] - x_ref)
    return J
end

c1 = DTO.Cost(obj1, nx, nu)
ct = DTO.Cost(objt, 2 * nx + 5, nu)
cT = DTO.Cost(objT, 2 * nx + 5, 0)
obj = [c1, [ct for t = 2:T-1]..., cT];

# ## constraints

ql = [-Inf; 0; -Inf; 0.0]
qu = [Inf; Inf; Inf; 1.0]
xl1 = [q1; ql]
xu1 = [q1; qu]
xlt = [ql; ql; -Inf * ones(5); -Inf * ones(nx)]
xut = [qu; qu; Inf * ones(5); Inf * ones(nx)]
ul = [-10.0; -10.0; zeros(nu - 2)]
uu = [10.0; 10.0; Inf * ones(nu - 2)]

bnd1 = DTO.Bound(nx, nu, xl=xl1, xu=xu1, ul=ul, uu=uu)
bndt = DTO.Bound(2 * nx + 5, nu, xl=xlt, xu=xut, ul=ul, uu=uu)
bndT = DTO.Bound(2 * nx + 5, 0, xl=xlt, xu=xut)
bnds = [bnd1, [bndt for t = 2:T-1]..., bndT];

function constraints_1(x, u, w)
    [
     # equality (8)
   	 RoboDojo.kinematics_foot(RoboDojo.hopper, x[1:RoboDojo.hopper.nq]) - RoboDojo.kinematics_foot(RoboDojo.hopper, x1[1:RoboDojo.hopper.nq]);
	 RoboDojo.kinematics_foot(RoboDojo.hopper, x[RoboDojo.hopper.nq .+ (1:RoboDojo.hopper.nq)]) - RoboDojo.kinematics_foot(RoboDojo.hopper, x1[RoboDojo.hopper.nq .+ (1:RoboDojo.hopper.nq)]);
     contact_constraints_equality(h, x, u, w);
     # inequality (12)
     contact_constraints_inequality_1(h, x, u, w);
    ]
end

function constraints_t(x, u, w)
    [
     # equality (4)
     contact_constraints_equality(h, x, u, w);
     # inequality (16)
     contact_constraints_inequality_t(h, x, u, w);
    ]
end

function constraints_T(x, u, w)
    x_travel = 0.5
    θ = x[nx + 5 .+ (1:nx)]
    [
     # equality (6)
	 x[1:RoboDojo.hopper.nq][collect([2, 3, 4])] - θ[1:RoboDojo.hopper.nq][collect([2, 3, 4])];
	 x[RoboDojo.hopper.nq .+ (1:RoboDojo.hopper.nq)][collect([2, 3, 4])] - θ[RoboDojo.hopper.nq .+ (1:RoboDojo.hopper.nq)][collect([2, 3, 4])];
     # inequality (10)
     x_travel - (x[1] - θ[1])
	 x_travel - (x[RoboDojo.hopper.nq + 1] - θ[RoboDojo.hopper.nq + 1])
     contact_constraints_inequality_T(h, x, u, w);
    ]
end

con1 = DTO.Constraint(constraints_1, nx, nu, idx_ineq=collect(8 .+ (1:12)))
cont = DTO.Constraint(constraints_t, 2nx + 5, nu, idx_ineq=collect(4 .+ (1:16)))
conT = DTO.Constraint(constraints_T, 2nx + 5, nu, idx_ineq=collect(6 .+ (1:10)))
cons = [con1, [cont for t = 2:T-1]..., conT];

# ## problem
p = DTO.solver(dyn, obj, cons, bnds,
    options=DTO.Options(
        tol=1.0e-2,
        constr_viol_tol=1.0e-2,
        print_level=0))

# ## initialize
x_interpolation = [x1, [[x1; zeros(5); x1] for t = 2:T]...]
u_guess = [[0.0; RoboDojo.hopper.gravity * RoboDojo.hopper.mass_body * 0.5 * h[1]; 1.0e-1 * ones(nu - 2)] for t = 1:T-1] # may need to run more than once to get good trajectory
DTO.initialize_states!(p, x_interpolation)
DTO.initialize_controls!(p, u_guess);

# ## solve
@time DTO.solve!(p);

# ## solution
x_sol, u_sol = DTO.get_trajectory(p)
@show x_sol[1]
@show x_sol[T]
sum([u[nu] for u in u_sol[1:end-1]])
x_sol[1] - x_sol[T][1:nx]

# ## visualize
vis = Visualizer()
render(vis)
q_sol = state_to_configuration([x[1:nx] for x in x_sol])
RoboDojo.visualize!(vis, RoboDojo.hopper, q_sol, Δt=h);

# ## comparison
function obj1_compare(x, u, w)
	x = x[1:8]
	u = u[1:2]

	J = 0.0
	J += 0.5 * transpose(x - x_ref) * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0]) * (x - x_ref)
	J += 0.5 * transpose(u) * Diagonal(r_cost * ones(hopper.nu)) * u
	return J
end

function objt_compare(x, u, w)
	x = x[1:8]
	u = u[1:2]
	J = 0.0
	J += 0.5 * transpose(x - x_ref) * q_cost * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0]) * (x - x_ref)
	J += 0.5 * transpose(u) * Diagonal(r_cost * ones(hopper.nu)) * u
	return J
end

function objT_compare(x, u, w)
	x = x[1:8]
	J = 0.0
	J += 0.5 * transpose(x - x_ref) * Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]) * (x - x_ref)
	return J
end

function obj_compare(x, u, w)
	J = 0.0
	J += obj1_compare(x[1], u[1], w[1])
	for t = 2:T-1
		J += objt_compare(x[t], u[t], w[t])
	end
	J += objT_compare(x[T], nothing, nothing)

	return J
end

obj_compare(x_sol, u_sol, [zeros(0) for t = 1:T])
@show norm(constraints_T(x_sol[T], zeros(0), zeros(0))[1:6], Inf)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Arrow, Tables, Test

include(joinpath(dirname(pathof(Arrow)), "../test/arrowjson.jl"))
# using .ArrowJSON

function runcommand(jsonname, arrowname, mode, verbose)
    if jsonname == ""
        error("must provide json file name")
    end
    if arrowname == ""
        error("must provide arrow file name")
    end

    if mode == "ARROW_TO_JSON"
        tbl = Arrow.Table(arrowname)
        df = ArrowJSON.DataFile(tbl)
        open(jsonname, "w") do io
            JSON3.write(io, df)
        end
    elseif mode == "JSON_TO_ARROW"
        df = ArrowJSON.parsefile(jsonname)
        open(arrowname, "w") do io
            Arrow.write(io, df)
        end
    elseif mode == "VALIDATE"
        df = ArrowJSON.parsefile(jsonname)
        tbl = Arrow.Table(arrowname)
        @test isequal(df, tbl)
    else
        error("unknown integration test mode: $mode")
    end
    return
end

max_forks_repo_name: XpressAI/arrow


text:
"""
    Keypoints(points, sz)
    Keypoints{N, T, M}(points, bounds)

`N`-dimensional keypoints represented as `SVector{N, T}`.

Spatial bounds are given by the polygon `bounds::Vector{SVector{N, T}}`
or `sz::NTuple{N, Int}`.

## Examples
{cell=Keypoints}
```julia
using DataAugmentation, StaticArrays
points = [SVector(y, x) for (y, x) in zip(4:5:80, 10:6:90)]
item = Keypoints(points, (100, 100))
```

{cell=Keypoints}
```julia
showitems(item)
```

"""
struct Keypoints{N, T, S<:Union{SVector{N, T}, Nothing}, M} <: AbstractArrayItem{M, S}
    data::AbstractArray{S, M}
    bounds::Bounds{N}
end

function Keypoints(data, sz::NTuple{N, Int}) where N
    return Keypoints(data, Bounds(sz))
end

Base.show(io::IO, item::Keypoints{N, T, M}) where {N, T, M} =
    print(io, "Keypoints{$N, $T, $M}() with $(length(item.data)) elements")

getbounds(keypoints::Keypoints) = keypoints.bounds

function project(P, keypoints::Keypoints{N, T}, bounds::Bounds{N}) where {N, T}
    # TODO: convert back to `T`?
    return Keypoints(
        map(fmap(P), keypoints.data),
        bounds,
    )
end
function project!(buf, P, keypoints::Keypoints{N, T}, bounds::Bounds{N}) where {N, T}
    # TODO: convert back to `T`?
    return Keypoints(
        map(fmap(P), keypoints.data),
        bounds,
    )
end

function showitem!(img, keypoints::Keypoints{N}) where N
    for point in filter(!isnothing, keypoints.data)
        showkeypoint!(img, point, RGBA(1, 0, 0, 1))
    end
    return img
end

# ## Wrappers
#
# We also define some wrappers for `Keypoints` which have the same representation,
# but a different meaning.

"""
    Polygon(points, sz)
    Polygon{N, T, M}(points, bounds)

Item wrapper around [`Keypoints`](#).

## Examples

{cell=Polygon}
```julia
using DataAugmentation, StaticArrays
points = [SVector(10., 10.), SVector(80., 20.), SVector(90., 70.), SVector(20., 90.)]
item = Polygon(points, (100, 100))
```

{cell=Polygon}
```julia
showitems(item)
```
"""
struct Polygon{N, T, S, M} <: ItemWrapper{Keypoints{N, T, S, M}}
    item::Keypoints{N, T, S, M}
end

Polygon(data, bounds) = Polygon(Keypoints(data, bounds))

Base.show(io::IO, item::Polygon{N, T, M}) where {N, T, M} =
    print(io, "Polygon{$N, $T}() with $(length(item.item.data)) elements")

function showitem!(img, polygon::Polygon)
    showpolygon!(img, itemdata(polygon), RGBA(1, 0, 0, 1))
    return img
end

"""
    BoundingBox(points, sz)
    BoundingBox{N, T, M}(points, bounds)

Item wrapper around [`Keypoints`](#).

## Examples

{cell=BoundingBox}
```julia
using DataAugmentation, StaticArrays
points = [SVector(10., 10.), SVector(80., 60.)]
item = BoundingBox(points, (100, 100))
```

{cell=BoundingBox}
```julia
showitems(item)
```
"""
struct BoundingBox{N, T, S} <: ItemWrapper{Keypoints{N, T, S, 1}}
    item::Keypoints{N, T, S, 1}
end

function BoundingBox(data::AbstractVector{<:SVector{N}}, bounds) where N
    length(data) == N || error("Give $N corner points for an $N-dimensional bounding box")
    BoundingBox(Keypoints(data, bounds))
end

Base.show(io::IO, item::BoundingBox{N, T}) where {N, T} =
    print(io, "BoundingBox{$N, $T}()")

function showitem!(img, bbox::BoundingBox{2})
    ul, br = itemdata(bbox)
    points = [
        ul,
        SVector(br[1], ul[2]),
        br,
        SVector(ul[1], br[2]),
    ]
    showpolygon!(img, points, RGBA(1, 0, 0, 1))
    return img
end

# ### Helpers

fmap(f) = x -> fmap(f, x)
fmap(f, x) = f(x)
fmap(f, x::Nothing) = x


text:
using Statistics

"""
    cv(estimator, X, y, splits; [verbosity = 1])

Cross-validate the `estimator` with features data `X` and label `y`. The iterable `splits` provides
vectors of indices for the training dataset. The remaining indices are used to create the
validation dataset.

Return a dictionary with an entry for the validation dataset and, if the parameter
`is_training_metric` is set in the `estimator`, an entry for the training dataset. Each entry of
the dictionary is another dictionary with an entry for each validation metric in the `estimator`.
Each of these entries is an array that holds the validation metric's value for each dataset, at the
last valid iteration.

# Arguments
* `estimator::LGBMEstimator`: the estimator to be fit.
* `X::Matrix{TX<:Real}`: the features data.
* `y::Vector{Ty<:Real}`: the labels.
* `splits`: the iterable providing arrays of indices for the training dataset.
* `verbosity::Integer`: keyword argument that controls LightGBM's verbosity. `< 0` for fatal logs
    only, `0` includes warning logs, `1` includes info logs, and `> 1` includes debug logs.
"""
function cv(estimator::LGBMEstimator, X::Matrix{TX}, y::Vector{Ty}, splits; verbosity::Integer = 1) where { TX<:Real,Ty<:Real }
    start_time = now()
    num_data = size(X)[1]
    ds_parameters = stringifyparams(estimator, DATASETPARAMS)
    bst_parameters = stringifyparams(estimator, BOOSTERPARAMS) #* " verbosity=$verbosity"

    full_ds = LGBM_DatasetCreateFromMat(X, ds_parameters)
    LGBM_DatasetSetField(full_ds, "label", y)

    split_scores = Dict{String,Dict{String,Vector{Float64}}}()
    for (split_idx, train_inds) in enumerate(splits)
        log_info(verbosity, "\nCross-validation: ", split_idx, "\n")

        log_debug(verbosity, "Started creating LGBM training dataset ", split_idx, "\n")
        train_ds = LGBM_DatasetGetSubset(full_ds, train_inds, ds_parameters)

        log_debug(verbosity, "Started creating LGBM test dataset ", split_idx, "\n")
        test_inds = setdiff(1:num_data, train_inds)
        test_ds = LGBM_DatasetGetSubset(full_ds, test_inds, ds_parameters)

        log_debug(verbosity, "Started creating LGBM booster ", split_idx, "\n")
        estimator.booster = LGBM_BoosterCreate(train_ds, bst_parameters)
        LGBM_BoosterAddValidData(estimator.booster, test_ds)

        results = train(estimator, ["validation"], verbosity, start_time)

        for dataset in keys(results)
            dataset_results = results[dataset]
            if !haskey(split_scores, dataset)
                split_scores[dataset] = Dict{String,Vector{Float64}}()
                for metric in keys(dataset_results)
                    split_scores[dataset][metric] = [dataset_results[metric][end]]
                end
            else
                for metric in keys(dataset_results)
                    push!(split_scores[dataset][metric], dataset_results[metric][end])
                end
            end
        end
    end

    log_info(verbosity, "\nCross-validation finished\n")
    cv_logsummary(split_scores, verbosity)

    return split_scores
end

function cv_logsummary(cv_results::Dict{String,Dict{String,Vector{Float64}}}, verbosity::Integer)
    for dataset in keys(cv_results)
        for metric in keys(cv_results[dataset])
            log_info(verbosity, "- ", dataset, "'s ", metric,
                     " mean: ", mean(cv_results[dataset][metric]),
                     ", std: ", std(cv_results[dataset][metric]), "\n")
        end
    end
end


text:
using RungeKutta

@testset "$(rpad("Tableau List",80))" begin

    tab_list = TableauList

    for tab in tab_list
        @test typeof(tab()) <: Tableau
    end

end


text:
module SystemTest

using Test
using Elsa, SparseArrays, LinearAlgebra
using Elsa: nlinks, nsites

@testset "operators" begin
    @test begin
        sys = System(LatticePresets.square(), Model(Hopping(@SMatrix[1 2; 2 1]))) |> grow(supercell = 3)
        h = sys.hamiltonian
        v = Elsa.boundaryoperator(h)
        nnz(h.matrix) == 144 && nnz(v.matrix) == 48 && v.inters == h.inters
    end
end

@testset "system" begin
    @test System(LatticePresets.honeycomb(), dim = Val(3), htype = Float64, ptype = Float32) isa System{3,2,Float32,Float64}
    @test nsites(System(LatticePresets.honeycomb()) |> grow(region = RegionPresets.square(300))) == 207946
    @test nsites(System(LatticePresets.square()) |> grow(supercell = 31)) == 961
    @test nsites(System(LatticePresets.bcc()) |> grow(region = RegionPresets.spheroid((10,4,4)))) == 1365

    @test nlinks(System(LatticePresets.honeycomb(), Model(Hopping(1, sublats = (1,2)))) |> grow(region = RegionPresets.square(30))) == 6094
    @test nlinks(System(LatticePresets.honeycomb(), Model(Hopping(1, ndists = (0,1)))) |> grow(region = RegionPresets.square(30))) == 3960
    @test nlinks(System(LatticePresets.honeycomb(), Model(Hopping(1, ndists = ((0,1),)))) |> grow(region = RegionPresets.square(30))) == 3960
    @test nlinks(System(LatticePresets.honeycomb(), Model(Hopping(1, sublats = (1,2), ndists = ((0,1), (0,1))))) |>
        grow(region = Region(:square, 30))) == 2040
    @test nlinks(System(LatticePresets.square(), Model(Hopping(1, range = 2))) |> grow(supercell = 31)) == 11532
    @test nlinks(System(LatticePresets.bcc(), Model(Hopping(1, range = 1))) |>
        grow(supercell = ((1, 2, 0),), region = RegionPresets.sphere(10))) == 10868
    @test nlinks(System(LatticePresets.bcc(), Model(Hopping(1, range = 1))) |>
        grow(supercell = (10, 20, 30), region = RegionPresets.sphere(4))) == 326
    @test nlinks(System(LatticePresets.bcc(), Model(Hopping(1, range = 1))) |> grow(supercell = (10, 20, 30))) == 84000

    @test nsites(SystemPresets.graphene_bilayer(twistindex = 31)) == 11908
end

@testset "system api" begin
    @test begin
        sys1 = System(LatticePresets.honeycomb(), Model(Hopping(1, sublats = (1,2))), dim = Val(3)) |> grow(region = RegionPresets.square(4)) |>
            transform(r -> r + SVector(0,0,1))
        sys2 = System(LatticePresets.square(), Model(Hopping(1, range = 2)), dim = Val(3)) |> grow(region = RegionPresets.square(4))
        sys3 = combine(sys1, sys2)
        (nlinks(sys3) == nlinks(sys1) + nlinks(sys2) == 294) &&
        (nsites(sys3) == nsites(sys1) + nsites(sys2) == 61)
    end

    @test begin
        sys1 = System(LatticePresets.honeycomb(), Model(Hopping(1, sublats = (1,2)))) |> grow(region = RegionPresets.circle(7))
        sys2 = System(LatticePresets.square(), Model(Hopping(1, range = 2))) |> grow(region = RegionPresets.circle(6))
        sys3 = combine(sys1, sys2, Model(Hopping(1, range = 1)))
        nlinks(sys3) == 5494
    end

    @test begin
        sys1 = System(LatticePresets.honeycomb())
        sys2 = transform(sys1, r -> 2r)
        sys2.lattice.sublats[1].sites[1] ≈ 2 * sys1.lattice.sublats[1].sites[1]
    end

    @test begin
        sys1 = System(LatticePresets.honeycomb())
        sys2 = System(LatticePresets.honeycomb())
        sys3 = System(LatticePresets.honeycomb())
        sys4 = combine(sys1, combine(sys2, sys3))
        allunique(sys4.sysinfo.names)
    end

    @test begin
        sys1 = System(LatticePresets.honeycomb(), Model(Onsite(@SMatrix[1 0; 0 2], sublats = :A),
                                        Onsite(@SMatrix[2 0; 0 3], sublats = :B)))
        sys2 = grow(sys1, region = RegionPresets.circle(1))
        sum(sys2.hamiltonian.matrix) == tr(sys2.hamiltonian.matrix) == 24
    end

    @test begin
        sys = System(LatticePresets.bcc(), Model(Hopping(1))) |> grow(supercell = 3) |> bound(except = 1)
        sys isa System{3,1}
    end

    @test begin
        sys = System(LatticePresets.cubic(), Model(Hopping(1, range = 3))) |> grow(supercell = (3,1,3)) |>
              bound(except = (1,3))
        sys isa System{3,2} && Elsa.nlinks(sys) == 252
    end

    @test begin
        sys = System(LatticePresets.honeycomb()) |> bound()
        sys isa System{2,0}
    end
end

# @test Elsa.nlinks(wrap(Lattice(:square, LinkRule(√2), Supercell(2)), exceptaxes = (1,))) == 14
# @test Elsa.nlinks(wrap(Lattice(:square, LinkRule(√2), Supercell(2)), exceptaxes = (2,))) == 14
# @test Elsa.nlinks(wrap(Lattice(:square, LinkRule(√2), Supercell(2)))) == 6

# @test begin
#     lat = mergesublats(Lattice(Preset(LatticePresets.honeycomb()_bilayer, twistindex = 2)), (2,1,1,1))
#     Elsa.nlinks(lat) == 32 && Elsa.nsublats(lat) == 2
# end

# @test begin
#     sys = System(Lattice(:square, LinkRule(2), Supercell(3)), Model(Onsite(1), Hopping(.3)))
#     nnz(sys.hbloch.matrix) == 81 && size(sys.vbloch.matrix) == (9, 9)
# end

# @test begin
#     sys = System(Lattice(LatticePresets.honeycomb(), LinkRule(1), Supercell(3)), Model(Onsite(1), Hopping(.3, (1,2))))
#     vel = velocity!(sys, k = (.2,.3))
#     ishermitian(vel) && size(vel) == (18, 18)
# end

# @test begin
#     sys = System(Lattice(LatticePresets.honeycomb(), LinkRule(1), Region(:square, 5)), Model(Onsite(1), Hopping(.3, (1,2))))
#     iszero(velocity!(sys))
# end

end # module


text:
#=
    Figure Snew
=#

using Inference
using Model

using JLD2
using Plots
using StatsPlots
using Random

include("../figure_defaults.jl")

##############################################################
## Load model and results
##############################################################

include("../../Results/MainResult_Setup.jl")
@load "Results/MainResult.jld2"

##############################################################
## Compute results
##############################################################

# Simulate model at 100 samples from the posterior
Θ = sample(C[1:100:end],100)
Y = [simulate_model(θ,T,param_dist,E;N=100_000,η) for θ in Θ]

##############################################################
## Figure 4(a,b)
##############################################################

#fig4ab = plot_fit_univariate(X,Y,T)
qlim=[0.0,3e4]
ulim=[0.0,150.0]

using KernelDensity

function vertical_density_multi!(plt,X1,X2,lim,i;fill=false,c=:blue)
    x = range(lim...,length=200)
    y1 = pdf(kde(X1),x)
    sc = maximum(y1)
    plot!(plt,i .+ 0.5 * y1/sc,x,fillrange=i,lc=c,fc=Model.brighten_colour(c,0.7),lw=1.0)
    for k = 1:length(X2)
        y2 = pdf(kde(X2[k]),x)
        plot!(plt,i .+ 0.5 * y2/sc,x,c="#333333",α=0.1,lw=1.0)
    end
end

plt1 = plot(legend=:none)
for i = 1:7
    vertical_density_multi!(plt1,X[i][2][1],[y[i][2][1] for y in Y],qlim,i+0.2,c=col_Q)
    vertical_density_multi!(plt1,X[i][1][1],[y[i][1][1] for y in Y],qlim,i,c=col_Q̄)
end
plot!(plt1,xticks=(1:7,T),box=:off,axis=:all,widen=true)

plt2 = plot(legend=:none)
for i = 1:7
    vertical_density_multi!(plt2,X[i][2][2],[y[i][2][2] for y in Y],ulim,i+0.2,c=col_U)
    vertical_density_multi!(plt2,X[i][1][2],[y[i][2][2] for y in Y],ulim,i,c=col_Ū)
end
plot!(plt2,xticks=(1:7,Int.(T)),box=:off,axis=:all,widen=true)

plot!(plt1,
    yformatter=y->y/1e4,
    ylabel="Cy5 (10⁴)",
    xticks=(1:7,[])
)

plot!(plt2,
    ylabel="BDP",
    xlabel="Time [min]"
)

plt = plot(plt1,plt2,layout=grid(2,1),size=(700,350))

plot!(plt,
    axis=:y,
    bottom_margin=2Plots.mm
)

savefig(plt,"$(@__DIR__)/figS5.svg")

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
type GlmResp{V<:FPVector} <: ModResp           # response in a glm model
    y::V                                       # response
    d::UnivariateDistribution
    l::Link
    devresid::V                                # (squared) deviance residuals
    eta::V                                     # linear predictor
    mu::V                                      # mean response
    mueta::V                                   # derivative of mu w.r.t. eta
    offset::V                                  # offset added to linear predictor (usually 0)
    var::V                                     # (unweighted) variance at current mu
    wts::V                                     # prior weights
    wrkwts::V                                  # working weights
    wrkresid::V                                # working residuals
    function GlmResp(y::V, d::UnivariateDistribution, l::Link,
                     eta::V, mu::V,
                     off::V, wts::V)
        if isa(d, Binomial)
            for yy in y; 0. <= yy <= 1. || error("$yy in y is not in [0,1]"); end
        else
            for yy in y; insupport(d, yy) || error("y must be in the support of d"); end
        end
        n = length(y)
        length(eta) == length(mu) == length(wts) == n || error("mismatched sizes")
        lo = length(off); lo == 0 || lo == n || error("offset must have length $n or length 0")
        res = new(y,d,l,similar(y),eta,mu,similar(y),off,similar(y),wts,similar(y),similar(y))
        updatemu!(res, eta)
        res
    end
end

# returns the sum of the squared deviance residuals
deviance(r::GlmResp) = sum(r.devresid)

# update the `devresid` field
devresid!(r::GlmResp) = devresid!(r.d, r.devresid, r.y, r.mu, r.wts)

# apply the link function generating the linear predictor (eta) vector from the mean vector (mu)
linkfun!(r::GlmResp) = linkfun!(r.l, r.eta, r.mu)

# apply the inverse link function generating the mean vector (mu) from the linear predictor (eta)
linkinv!(r::GlmResp) = linkinv!(r.l, r.mu, r.eta)

# evaluate the mueta vector (derivative of mu w.r.t. eta) from the linear predictor (eta)
mueta!(r::GlmResp) = mueta!(r.l, r.mueta, r.eta)

function updatemu!{T<:FPVector}(r::GlmResp{T}, linPr::T)
    n = length(linPr)
    length(r.offset) == n ? map!(Add(), r.eta, linPr, r.offset) : copy!(r.eta, linPr)
    linkinv!(r); mueta!(r); var!(r); wrkresid!(r); devresid!(r)
    sum(r.devresid)
end

updatemu!{T<:FPVector}(r::GlmResp{T}, linPr) = updatemu!(r, convert(T,vec(linPr)))

var!(r::GlmResp) = var!(r.d, r.var, r.mu)

wrkresid!(r::GlmResp) = map1!(Divide(), map!(Subtract(), r.wrkresid, r.y, r.mu), r.mueta)

function wrkresp(r::GlmResp)
    if length(r.offset) > 0 return map1!(Add(), map(Subtract(), r.eta, r.offset), r.wrkresid) end
    map(Add(), r.eta, r.wrkresid)
end

function wrkwt!(r::GlmResp)
    wrkwts = r.wrkwts
    mueta = r.mueta
    var = r.var
    if length(r.wts) == 0
        for i = 1:length(r.var)
            @inbounds wrkwts[i] = abs2(mueta[i])/var[i]
        end
    else
        wts = r.wts
        for i = 1:length(r.var)
            @inbounds wrkwts[i] = wts[i] * abs2(mueta[i])/var[i]
        end
    end
    wrkwts
end

type GeneralizedLinearModel <: LinPredModel
    rr::GlmResp
    pp::LinPred
    fit::Bool
end

function coeftable(mm::GeneralizedLinearModel)
    cc = coef(mm)
    se = stderr(mm)
    zz = cc ./ se
    CoefTable(hcat(cc,se,zz,2.0 * ccdf(Normal(), abs(zz))),
              ["Estimate","Std.Error","z value", "Pr(>|z|)"],
              ["x$i" for i = 1:size(mm.pp.X, 2)], 4)
end

function confint(obj::GeneralizedLinearModel, level::Real)
    hcat(coef(obj),coef(obj)) + stderr(obj)*quantile(Normal(),(1. -level)/2.)*[1. -1.]
end
confint(obj::GeneralizedLinearModel) = confint(obj, 0.95)

deviance(m::GeneralizedLinearModel)  = deviance(m.rr)

function StatsBase.fit(m::GeneralizedLinearModel; verbose::Bool=false, maxIter::Integer=30,
                       minStepFac::Real=0.001, convTol::Real=1.e-6)
    m.fit && return m
    maxIter >= 1 || error("maxIter must be positive")
    0 < minStepFac < 1 || error("minStepFac must be in (0, 1)")

    cvg = false; p = m.pp; r = m.rr
    scratch = similar(p.X)
    devold = updatemu!(r, linpred(delbeta!(p, wrkresp(r), wrkwt!(r), scratch)))
    installbeta!(p)
    for i=1:maxIter
        f = 1.0
        dev = updatemu!(r, linpred(delbeta!(p, r.wrkresid, wrkwt!(r), scratch)))
        while dev > devold
            f /= 2.; f > minStepFac || error("step-halving failed at beta0 = $(p.beta0)")
            dev = updatemu!(r, linpred(p, f))
        end
        installbeta!(p, f)
        crit = (devold - dev)/dev
        verbose && println("$i: $dev, $crit")
        if crit < convTol; cvg = true; break end
        devold = dev
    end
    cvg || error("failure to converge in $maxIter iterations")
    m.fit = true
    m
end

function StatsBase.fit(m::GeneralizedLinearModel, y; wts=nothing, offset=nothing, dofit::Bool=true, fitargs...)
    r = m.rr
    V = typeof(r.y)
    r.y = isa(y, V) ? copy(y) : convert(V, y)
    wts == nothing || (r.wts = isa(wts, V) ? copy(wts) : convert(V, wts))
    offset == nothing || (r.offset = isa(offset, V) ? copy(offset) : convert(V, offset))
    r.mu = mustart(r.d, r.y, r.wts)
    linkfun!(r.l, r.eta, r.mu)
    updatemu!(r, r.eta)
    fill!(m.pp.beta0, zero(eltype(m.pp.beta0)))
    m.fit = false
    dofit ? fit(m; fitargs...) : m
end

function StatsBase.fit{T<:FloatingPoint,V<:FPVector}(::Type{GeneralizedLinearModel},
                                                     X::Matrix{T}, y::V, d::UnivariateDistribution,
                                                     l::Link=canonicallink(d);
                                                     dofit::Bool=true,
                                                     wts::V=fill!(similar(y), one(eltype(y))),
                                                     offset::V=similar(y, 0), fitargs...)
    size(X, 1) == size(y, 1) || DimensionMismatch("number of rows in X and y must match")
    n = length(y)
    length(wts) == n || error("length(wts) = $lw should be 0 or $n")
    wts = T <: Float64 ? copy(wts) : convert(typeof(y), wts)
    off = T <: Float64 ? copy(offset) : convert(Vector{T}, offset)
    mu = mustart(d, y, wts)
    eta = linkfun!(l, similar(mu), mu)
    if !isempty(off)
        subtract!(eta, off)
    end
    rr = GlmResp{typeof(y)}(y, d, l, eta, mu, offset, wts)
    res = GeneralizedLinearModel(rr, DensePredChol(X), false)
    dofit ? fit(res; fitargs...) : res
end

StatsBase.fit(::Type{GeneralizedLinearModel}, X::Matrix, y::AbstractVector, d::UnivariateDistribution,
              l::Link=canonicallink(d); kwargs...) =
    fit(GeneralizedLinearModel, float(X), float(y), d, l; kwargs...)

glm(X, y, args...; kwargs...) = fit(GeneralizedLinearModel, X, y, args...; kwargs...)

## scale(m) -> estimate, s, of the scale parameter
## scale(m,true) -> estimate, s^2, of the squared scale parameter
type DispersionFun <: Functor{2} end
evaluate{T<:FP}(::DispersionFun,wt::T,resid::T) = wt*abs2(resid)
function scale(m::GeneralizedLinearModel, sqr::Bool=false)
    if isa(m.rr.d, Union(Binomial, Poisson))
        return 1.
    end

    s = sum(DispersionFun(), m.rr.wrkwts, m.rr.wrkresid)/df_residual(m)
    sqr ? s : sqrt(s)
end


text:
# ---
# title: 315. Count of Smaller Numbers After Self
# id: problem315
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search, Divide and Conquer, Sort, Binary Indexed Tree, Segment Tree
# link: <https://leetcode.com/problems/count-of-smaller-numbers-after-self/description/>
# hidden: true
# ---
#
# You are given an integer array _nums_ and you have to return a new _counts_
# array. The _counts_ array has the property where `counts[i]` is the number of
# smaller elements to the right of `nums[i]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [5,2,6,1]
#     Output: [2,1,1,0]
#     Explanation:
#     To the right of 5 there are **2** smaller elements (2 and 1).
#     To the right of 2 there is only **1** smaller element (1).
#     To the right of 6 there is **1** smaller element (1).
#     To the right of 1 there is **0** smaller element.
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 10^5`
#   * `-10^4 <= nums[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

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


text:
import Pkg
Pkg.activate("/code/")
using KeboolaConnectionTransformation

KeboolaConnectionTransformation.run()

max_stars_repo_name: keboola/julia-transformation

max_issues_repo_name: keboola/julia-transformation

max_forks_repo_name: keboola/julia-transformation


text:
import AlgebraicMultigrid: Level, MultiLevel, GaussSeidel

function multigrid(A::TA; max_levels = 10, max_coarse = 10,
                    presmoother = GaussSeidel(), postsmoother = GaussSeidel()) where {T,V,TA<:SparseMatrixCSC{T,V}}

    levels = Vector{Level{TA,TA,TA}}()
    w = AlgebraicMultigrid.MultiLevelWorkspace(Val{1}, eltype(A))

    while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
        AlgebraicMultigrid.residual!(w, size(A, 1))
        A = extend!(levels, A)
        AlgebraicMultigrid.coarse_x!(w, size(A, 1))
        AlgebraicMultigrid.coarse_b!(w, size(A, 1))
        #=if size(A, 1) <= max_coarse
            # push!(levels, Level(A,spzeros(T,0,0),spzeros(T,0,0)))
            break
        end=#
    end

    MultiLevel(levels, A, Pinv(A), presmoother, postsmoother, w)
end

function extend!(levels, A::SparseMatrixCSC{Ti,Tv}) where {Ti,Tv}

    size_F = size(A, 1)
    size_C = rem(size_F,2) == 0 ? div((size_F-1), 2) +1 : div((size_F-1), 2)
    total_len = size_C + 2 * (size_C - 1)

    I = Vector{Tv}(undef, total_len)
    J = Vector{Tv}(undef, total_len)
    V = Vector{Ti}(undef, total_len)

    l = 1
    for k = 1:size_C
        I[l], J[l], V[l] = 2k, k, 1
        l += 1
    end
    for k = 1:size_C - 1
        I[l], J[l], V[l] = 2k+1, k, 0.5
        I[l+1], J[l+1], V[l+1] = 2k+1, k+1, 0.5
        l += 2
    end

    P = sparse(I, J, V, size_F, size_C)

    R = AlgebraicMultigrid.adjoint(P)

    push!(levels, Level(A, P, R))

    R * A * P
end

Base.show(io::IO, level::Level) = print(io, "Level")


text:
"""
    DoubleFloats

A package for extended precision floating point math.

exported types: `Double64`, `Double32`, `Double16`
                `ComplexDF64`, `ComplexDF32`, `ComplexDF16`
"""
module DoubleFloats

abstract type MultipartFloat <: AbstractFloat end

export DoubleFloat,
       Double64, Double32, Double16,
       ComplexDF64, ComplexDF32, ComplexDF16,
       @df64_str, @df32_str, @df16_str,
       @via_df64, @via_df32,
       MultipartFloat, HI, LO, HILO,
       stringtyped, showtyped, showall,
       isnonzero, ispositive, isnegative, isnonpositive, isnonnegative,
       isposinf, isneginf, isfractional, isnormal, issquare,
       nan, inf, posinf, neginf,
       ulp,
       intpart, fracpart, fmod,
       square, cube,
       add2, sub2, mul2, div2,
       ⊕, ⊖, ⊗, ⊘,
       spread, sld, tld,
       signs,
       randpm,
       tanpi,
       agm, agm1, ellipk,
       besselj0, besselj1, besselj, bessely0, bessely1, bessely,
       gamma, loggamma, lgamma, erf, erfc,
       matrixfunction

using Base.MathConstants: pi, golden, ℯ, eulergamma, catalan

using Requires
using Printf
using Random
import Random: rand, randn

using Polynomials

using GenericSVD, GenericSchur
using LinearAlgebra
import LinearAlgebra: norm, mul!

import Quadmath

import Base: hash, promote_type, promote_rule, convert,
             string, show, parse, tryparse, eltype,
             signbit, sign, abs, flipsign, copysign,
             significand, exponent, precision, widen,
             (+), (-), (*), (/), (\), (^), inv, sqrt, cbrt,
             (==), (!=), (<), (<=), (>=), (>), isequal, isless,
             IEEEFloat, iszero, isone, isinf, isnan, isinf, isfinite, issubnormal,
             isinteger, isodd, iseven, zero, one,
             typemax, typemin, floatmax, floatmin, maxintfloat,
             min, max, minmax, minimum, maximum,
             floor, ceil, trunc, round, div, fld, cld,
             rem, mod, rem2pi, mod2pi, divrem, fldmod,
             sqrt, cbrt,
             BigFloat, BigInt,
             Int8, Int16, Int32, Int64, Int128,
             Float64, Float32, Float16

import Base.Math: modf, fma, muladd,
             log, log1p, log2, log10, exp, expm1, exp2, exp10,
             sin, cos, tan, csc, sec, cot, sincos, sinpi, cospi,
             asin, acos, atan, acsc, asec, acot,
             sinh, cosh, tanh, csch, sech, coth,
             asinh, acosh, atanh, acsch, asech, acoth

import Quadmath: Float128

include("Double.jl")   # Double64, Double32, Double16

include("math/errorfree.jl")
include("math/leasterror.jl")

include("doubletriple/double.jl")
include("doubletriple/double_consts.jl")
include("doubletriple/triple.jl")
include("doubletriple/triple_consts.jl")
include("doubletriple/triple_pi.jl")

include("type/constructors.jl")
include("type/promote.jl")
include("type/convert.jl")
include("type/compare.jl")
include("type/predicates.jl")
include("type/specialvalues.jl")
include("type/string.jl")
include("type/show.jl")
include("type/parse.jl")
include("type/quadmath.jl")

include("math/prearith/prearith.jl")
include("math/prearith/floorceiltrunc.jl")
include("math/prearith/divrem.jl")

include("math/arithmetic/fma.jl")
include("math/arithmetic/modpi.jl")
include("math/arithmetic/normalize_hypot.jl")
include("math/arithmetic/mixedarith.jl")
include("math/arithmetic/complex/pow.jl")
include("math/ops.jl")

include("math/elementary/sequences.jl")
include("math/elementary/explog.jl")
include("math/elementary/trig.jl")
include("math/elementary/arctrig.jl")
include("math/elementary/hyptrig.jl")
include("math/elementary/archyptrig.jl")
include("math/elementary/templated.jl")
include("math/elementary/complex.jl")

include("math/linearalgebra/matmul.jl")
include("math/linearalgebra/support.jl")
include("math/linearalgebra/matrixfunction.jl")

include("extras/random.jl")
include("extras/misc.jl")
include("extras/macros.jl")

if VERSION >= v"1.1.0"
  using SpecialFunctions
  include("math/special/specialfunctions.jl")
else
    function __init__()
        @require SpecialFunctions="276daf66-3868-5448-9aa4-cd146d93841b" begin
            include("math/special/specialfunctions.jl")
        end
    end
end

end # module DoubleFloats


text:
using Test

include("luhn.jl")

@testset "invalid" begin
    @testset "single digit strings can not be valid" begin
        @test !luhn("1")
    end

    @testset "A single zero is invalid" begin
        @test !luhn("0")
    end

    @testset "invalid Canadian SIN" begin
        @test !luhn("046 454 287")
    end

    @testset "invalid credit card" begin
        @test !luhn("8273 1232 7352 0569")
    end

    @testset "valid strings with a non-digit added become invalid" begin
        @test !luhn("046a 454 286")
    end
end

@testset "valid" begin
    @testset "valid Canadian SIN" begin
        @test luhn("046 454 286")
    end
end


text:
module GRF

# using ARCHModels
using Cbc
# using CSV
using Dates
using DataFrames
using Distributions
# using HypothesisTests
using HTTP
using Ipopt
using JSON
using JuMP
using Juniper
# using LinearAlgebra
# using Random
using Statistics
using StatsPlots

include("modulo01.jl")

# using Pkg
# println(haskey(Pkg.installed(), "JuMP"))

#=
nikkein = CSV.read("NIKKEI225n.csv", DataFrame)
println(nikkein[!,:Simbolo])
=#

#=
nikkeic = yahoo(NIKKEI225n, Date(2015,12,31), Date(2021,3,31))
println(describe(nikkeic))
CSV.write("NIKKEI225c.csv", nikkeic)
nikkeir = retornos(nikkeic, :log)
CSV.write("NIKKEI225r.csv", nikkeir)
=#

#=
mibn = CSV.read("FTSEMIBn.csv", DataFrame)
println(mibn[!,:Simbolo])
=#

#=
mibc = yahoo(FTSEMIBn, Date(2015,12,31), Date(2021,3,31))
println(describe(mibc))
CSV.write("FTSEMIBc.csv", mibc)
mibr = retornos(mibc, :log)
CSV.write("FTSEMIBr.csv", mibr)
=#

# Notebooks interativos
# binder: https://mybinder.org
# https://mybinder.org/v2/gh/ASaragga/GRF.jl/HEAD
#
# nbviewer: https://nbviewer.jupyter.org
# https://github.com/ASaragga/GRF.jl/blob/master/src/Modulo1.1-1.7.ipynb

end # module


text:
@graphs begin
	def_card_prio = """
	After the Free Peoples' player have selected their combat card, select and play a card. If no card is matching the 4 first items in the priority list, do not play any card.

	Priority:
	1. Strategy card that cancels the Free Peoples' card
	2. Doesn't use the term "Fellowship revealed"
	3. Doesn't add a hunt tile or corruption
	4. Character card
	5. Ascending order of initiative
	6. Random

	If the selected card requires units to be downgraded or sacrificed, do so as long as the army does not become non-*aggressive*.
	"""

	sortie_card_prio = """
	After the Free Peoples' player have selected their combat card, select and play a card. If no card is matching the 2 first items in the priority list, do not play any card.

	Priority:
	1. Charcter card that doesn't use the term "Fellowship revealed"
	2. Character card that doesn't add a hunt tile or corruption
	3. Ascending order of initiative
	4. Random

	If the selected card requires units to be downgraded or sacrificed, do so as long as the army does not become non-*aggressive*.
	"""

	wk_card_prio = """
	After the Free Peoples' player have selected their combat card, select and play a card.

	Priority:
	1. Strategy card
	2. Durin's Bane
	3. Character card
	4. Doesn't use the term "Fellowship revealed"
	5. Doesn't add a hunt tile or corruption
	6. Ascending order of initiative
	7. Random

	If the selected card requires units to be downgraded or sacrificed, do so as long as the army does not become non-*aggressive*.
	"""

	attack_card_prio = """
	After the Free Peoples' player have selected their combat card, select and play a card.

	Priority:
	1. Durin's Bane
	2. Strategy card
	3. Character card
	4. Doesn't use the term "Fellowship revealed"
	5. Doesn't add a hunt tile or corruption
	6. Ascending order of initiative
	7. Random

	If the selected card requires units to be downgraded or sacrificed, do so as long as the army does not become non-*aggressive*.
	"""

	################################################################################
	 @node battle = Start() -> rearguard
	 @node rearguard = PerformAction("All units from nations not at war form the rearguard.") -> army_attacking
	 @node army_attacking = BinaryCondition("The Shadow army is attacking.") -> [n_true = is_sortie, n_false = def_in_stronghold]

	 ########################################
	 @node def_in_stronghold = BinaryCondition("""
					 The Shadow army is defending in a region with a stronghold.
					 """) -> [n_true = should_retreat_to_stronghold, n_false = field_def_card_prio]
	 @node field_def_card_prio = PerformAction(def_card_prio) -> field_def_resolve
	 @node field_def_resolve = JumpToGraph("battle_resolve") -> field_attacking_fp_continues
	 @node field_attacking_fp_continues = BinaryCondition("""
														  The Free Peoples' player is continuing the attack.
														  """) -> [n_true = retreat_prio, n_false = field_def_end]
	 @node retreat_prio = PerformAction("""
				   Retreat from combat into region accroding to the following prioirty.

				   Priority:
				   1. Does not creat a *threat*
				   2. Reduce distance to *target* or *exposed* region
				   3. Increase the number of *mobile* armies
				   4. Increase the number of *aggressive* armies
				   5. Contains a settlement
				   6. Contains the highest *value* army
				   7. Adjacent to the highest *value* army
				   8. Random
				   """) -> field_def_end
	 @node field_def_end = End("End of Battle") -> []

	 ########################################
	 @node should_retreat_to_stronghold = BinaryCondition("""
					 The Shadow army is not under siege.
					 And, the *value* is less or equal to the attacking army's.
					 And, the number of units is less than 8.
					 """) -> [n_true = retreat_to_stronghold, n_false = def_card_prio]
	 @node def_card_prio = PerformAction(def_card_prio) -> def_resolve
	 @node def_resolve = JumpToGraph("battle_resolve") -> attacking_fp_continues
	 @node attacking_fp_continues = BinaryCondition("""
													The Free Peoples' player is continuing the attack.
													""") -> [n_true = should_retreat_to_stronghold, n_false = def_end]
	 @node def_end = End("End of Battle") -> []

	 @node retreat_to_stronghold = PerformAction("Retreat into stronghold.") -> retreat_stronghold_end
	 @node retreat_stronghold_end = End("End of Battle") -> []

	 ########################################
	 @node is_sortie = BinaryCondition("Battle is a sortie.") -> [n_true = sortie_card_prio, n_false = army_with_wk]
	 @node sortie_card_prio = PerformAction(sortie_card_prio) -> sortie_resolve
	 @node sortie_resolve = JumpToGraph("battle_resolve") -> sortie_round_end
	 @node sortie_round_end = JumpToGraph("battle_round_end") -> sortie_card_prio

	 ########################################
	 @node army_with_wk = BinaryCondition("Army include the Witch King.") -> [n_true = wk_card_prio, n_false = should_play_card]
	 @node wk_card_prio = PerformAction(wk_card_prio) -> wk_resolve
	 @node wk_resolve = JumpToGraph("battle_resolve") -> wk_round_end
	 @node wk_round_end = JumpToGraph("battle_round_end") -> should_play_card

	 ########################################
	 @node should_play_card = BinaryCondition("""
					 The Shadow is conducting a siege.
					 Or, the Shadow is holding more than 4 cards.
					 """) -> [n_true = attack_card_prio, n_false = attack_play_no_card]

	 @node attack_card_prio = PerformAction(attack_card_prio) -> attack_resolve
	 @node attack_play_no_card = PerformAction("Do not play a combat card.") -> attack_resolve
	 @node attack_resolve = JumpToGraph("battle_resolve") -> attack_round_end
	 @node attack_round_end = JumpToGraph("battle_round_end") -> should_play_card

	 ################################################################################
	 @node battle_resolve = Start() -> roll
	 @node roll = PerformAction("Roll for combat and reroll misses.") -> casualties
	 @node casualties = PerformAction("""
				   Remove casualties.

				   Priority:
				   1. Maximizes effect of the card played
				   2. Retains highest army *value* with lowest number of units
				   3. Keeps one unit of each nation
				   4. Random
				   """) -> battle_resolve_return
	 @node battle_resolve_return = ReturnFromGraph() -> []

	 ################################################################################
	 @node battle_round_end = Start() -> is_fp_dead
	 @node is_fp_dead = BinaryCondition("There are Free Peoples' units remaining.") -> [n_true = press_on, n_false = no_fp_left]
	 @node no_fp_left = BinaryCondition("""
					 Moving in to the conquered region would:
					 - win the game; or
					 - decrease distance to *target*; or
					 - remove a *threat*.
					 """) -> [n_true = move_into_conquered, n_false = end_without_moving]

	 @node move_into_conquered = PerformAction("Move the largest *value* possible into the conquered region.") -> move_into_conquered_end
	 @node move_into_conquered_end = End("End of Battle") -> []

	 @node end_without_moving = PerformAction("Do not move any units into the conquered region.") -> end_without_moving_end
	 @node end_without_moving_end = End("End of Battle") -> []

	 @node press_on = BinaryCondition("""
					 A field battle was fought.
					 """) -> [n_true = aggressive_if_continue, n_false = mili_strat]
	 @node mili_strat = CheckStrategy("military") -> [n_true = aggressive_if_continue, n_false = press_on_2]
	 @node press_on_2 = BinaryCondition("""
					 The Fellowship is on the Mordor track.
					 """) -> [n_true = another_round_if_possible, n_false = no_more_round]
	 @node aggressive_if_continue = BinaryCondition("""
					 The Shadow army is *aggressive* and, if a siege battle is fought, would remain *aggressive* after an Elite downgrade to continue the battle.
					 """) -> [n_true = another_round_if_possible, n_false = no_more_round_2]
	 @node another_round_if_possible = BinaryCondition("A siege battle is fought and the Shadow army have no Elites left") -> [n_true = no_more_round_2, n_false = one_more_round]
	 @node one_more_round = PerformAction("Continue the battle, downgrade an Elite if necessary.") -> one_more_round_return
	 @node one_more_round_return = ReturnFromGraph() -> []

	 @node no_more_round = PerformAction("End battle") -> no_more_round_end
	 @node no_more_round_end = End("End of Battle") -> []

	 @node no_more_round_2 = PerformAction("End battle") -> no_more_round_end_2
	 @node no_more_round_end_2 = End("End of Battle") -> []

end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_CONTEXT, SA_REORDER, SA_REPLACE_CALLS)

function pagerank(A, p, r, d_inv, maxiter) # p: initial rank, r: damping factor
#  set_matrix_property(Dict(
#    :A => SA_SYMM_STRUCTURED | SA_SYMM_VALUED | SA_STRUCTURE_ONLY))

  bytes = maxiter*(nnz(A)*4 + size(A, 1)*4*8)

  p = p.*d_inv
  d_inv = copy(d_inv)
  Ap = zeros(size(A, 1))

  t = time()

  for i = 1:maxiter
    Ap = (1-r)*A*p + r
    Ap = Ap.*d_inv

    if i == maxiter
      err = norm(Ap - p)/norm(p)
      println("error = $err")
    end

    temp = Ap
    Ap = p
    p = temp
  end

  t = time() - t
  println("pagerank takes $t sec ($(bytes/t/1e9) gbps)")

  p
end

function pagerank_with_init(A0, r, d_inv, maxiter) # p: initial rank, r: damping factor
  set_matrix_property(Dict(
    :A0 => SA_SYMM_VALUED)
  )

  A0 = spones(A0)

  m = size(A0, 1)
  p = repmat([1/m], m)

  d = max(convert(Array{eltype(A0),1}, vec(sum(A0, 2))), 1) # num of neighbors
  A = scale(A0,1./d)

  bytes = maxiter*(nnz(A)*4 + m*4*8)

  p = p.*d_inv
  d_inv = copy(d_inv)
  Ap = zeros(size(A, 1))

  t = time()

  for i = 1:maxiter
    Ap = (1-r)*A*p + r
    Ap = Ap.*d_inv

    if i == maxiter
      err = norm(Ap - p)/norm(p)
      println("error = $err")
    end

    temp = Ap
    Ap = p
    p = temp
  end

  t = time() - t
  println("pagerank takes $t sec ($(bytes/t/1e9) gbps)")

  p
end

function pagerank_call_replacement(A, p, r, d_inv, maxiter, do_print) # p: initial rank, r: damping factor
  bytes = maxiter*(nnz(A)*4 + size(A, 1)*4*8)
  A = copy(A)
  p = p.*d_inv
  d_inv = copy(d_inv)
  Ap = copy(p)

  t = time()
  reorder_time = 0.

  __mknobA = (Sparso.new_matrix_knob)(:A, true, true, true, true, false, false)

  fknob_spmv = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknobA, fknob_spmv)

  for i = 1:maxiter
    #Ap = ((1-r) *A * p + r).*d_inv
    Sparso.SpMV!(Ap, 1 - r, A, p, 0, p, r, d_inv, fknob_spmv)
    if do_print && i == maxiter
      err = norm(Ap - p)/norm(p)
      println("error = $err")
    end

    temp = Ap
    Ap = p
    p = temp
  end

  (Sparso.delete_matrix_knob)(__mknobA)
  (Sparso.delete_function_knob)(fknob_spmv)

  t = time() - t
  if do_print
    println("pagerank takes $t sec ($(bytes/t/1e9) gbps)")
  end

  p
end

function pagerank_call_replacement_and_context_opt(A, p, r, d_inv, maxiter, do_print) # p: initial rank, r: damping factor
  bytes = maxiter*(nnz(A)*4 + size(A, 1)*4*8)
  A = copy(A)
  p = p.*d_inv
  d_inv = copy(d_inv)
  Ap = copy(p)

  t = time()
  reorder_time = 0.

  __mknobA = (Sparso.new_matrix_knob)(:A, true, true, true, true, true, false)

  fknob_spmv = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknobA, fknob_spmv)

  for i = 1:maxiter
    #Ap = ((1-r) *A * p + r).*d_inv
    Sparso.SpMV!(Ap, 1 - r, A, p, 0, p, r, d_inv, fknob_spmv)
    if do_print && i == maxiter
      err = norm(Ap - p)/norm(p)
      println("error = $err")
    end

    temp = Ap
    Ap = p
    p = temp
  end

  (Sparso.delete_matrix_knob)(__mknobA)
  (Sparso.delete_function_knob)(fknob_spmv)

  t = time() - t
  if do_print
    println("pagerank takes $t sec ($(bytes/t/1e9) gbps)")
  end

  p
end

function pagerank_call_replacement_and_context_opt_and_reordering(A, p, r, d_inv, maxiter, do_print) # p: initial rank, r: damping factor
  bytes = maxiter*(nnz(A)*4 + size(A, 1)*4*8)
  A = copy(A)
  p = p.*d_inv
  d_inv = copy(d_inv)
  Ap = copy(p)

  t = time()
  reorder_time = 0.

  __mknobA = (Sparso.new_matrix_knob)(:A, true, true, true, true, true, false)

  fknob_spmv = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknobA, fknob_spmv)

  (Sparso.set_reordering_decision_maker)(fknob_spmv)
  reordering_status = [false, C_NULL, C_NULL, C_NULL, C_NULL, reorder_time]

  for i = 1:maxiter
    #Ap = ((1-r) *A * p + r).*d_inv
    Sparso.SpMV!(Ap, 1 - r, A, p, 0, p, r, d_inv, fknob_spmv)
    if do_print && i == maxiter
      err = norm(Ap - p)/norm(p)
      println("error = $err")
    end

    temp = Ap
    Ap = p
    p = temp
  end

  (Sparso.delete_matrix_knob)(__mknobA)
  (Sparso.delete_function_knob)(fknob_spmv)
  (Sparso.reverse_reordering)(reordering_status, :__delimitor__, p, Sparso.ROW_PERM)

  t = time() - t
  if do_print
    println("pagerank takes $t sec ($(bytes/t/1e9) gbps)")
  end

  p
end

A0 = matrix_market_read(ARGS[1], true, true)
A0 = spones(A0)

m = size(A0, 1)
p = repmat([1/m], m)
r = 0.15

d = max(convert(Array{eltype(A0),1}, vec(sum(A0, 2))), 1) # num of neighbors
A = scale(A0,1./d)

maxiter = 100
bytes = maxiter*(nnz(A)*4 + m*4*8)

d_inv = 1./d
x = pagerank(A0, p, r, d_inv, maxiter)
println("\nOriginal: ")
x = pagerank(A0, p, r, d_inv, maxiter)
println("End original.")

println("\nManual_call_replacement:")
x = pagerank_call_replacement(A0, p, r, d_inv, maxiter, false)

Sparso.reset_spmp_spmv_time()
Sparso.reset_knob_spmv_time()
x = pagerank_call_replacement(A0, p, r, d_inv, maxiter, true)
t = Sparso.get_spmp_spmv_time()
println("time spent on spmp spmv $t sec ($(bytes/t/1e9) gbps)")
t = Sparso.get_knob_spmv_time()
println("time spent on knob spmv $t sec ($(bytes/t/1e9) gbps)")
println("End Manual_call_replacement.")

println("\nManual_call_replacement_and_context_opt:")
x = pagerank_call_replacement_and_context_opt(A, p, r, d_inv, maxiter, false)

Sparso.reset_spmp_spmv_time()
Sparso.reset_knob_spmv_time()
#Sparso.set_knob_log_level(1)
x = pagerank_call_replacement_and_context_opt(A0, p, r, d_inv, maxiter, true)
t = Sparso.get_spmp_spmv_time()
println("time spent on spmp spmv $t sec ($(bytes/t/1e9) gbps)")
t = Sparso.get_knob_spmv_time()
println("time spent on knob spmv $t sec ($(bytes/t/1e9) gbps)")

println("End Manual_call_replacement_and_context_opt.")

println("\nManual_call_replacement_and_context_opt_and_reorder:")
x = pagerank_call_replacement_and_context_opt_and_reordering(A, p, r, d_inv, maxiter, false)

Sparso.reset_spmp_spmv_time()
Sparso.reset_knob_spmv_time()
#Sparso.set_knob_log_level(1)
x = pagerank_call_replacement_and_context_opt_and_reordering(A0, p, r, d_inv, maxiter, true)
t = Sparso.get_spmp_spmv_time()
println("time spent on spmp spmv $t sec ($(bytes/t/1e9) gbps)")
t = Sparso.get_knob_spmv_time()
println("time spent on knob spmv $t sec ($(bytes/t/1e9) gbps)")

println("End Manual_call_replacement_and_context_opt_and_reorder.")

# copy A since we change A in-place
#Sparso.set_knob_log_level(0)
A2 = copy(A0)
#@acc x= pagerank(A2, p, r, d_inv, maxiter)
@acc x= pagerank_with_init(A2, r, d_inv, maxiter)

println("\nAccelerated: ")
Sparso.reset_spmp_spmv_time()
Sparso.reset_knob_spmv_time()
#Sparso.set_knob_log_level(1)
# It seems that p has been changed by sparse accelerator
p = repmat([1/m], m)
#@acc x = pagerank(A0, p, r, d_inv, maxiter)
@acc x = pagerank_with_init(A2, r, d_inv, maxiter)
t = Sparso.get_spmp_spmv_time()
println("time spent on spmp spmv $t sec ($(bytes/t/1e9) gbps)")
t = Sparso.get_knob_spmv_time()
println("time spent on knob spmv $t sec ($(bytes/t/1e9) gbps)")
println("End accelerated.")

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso