text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

import ..SimData
import ..rawdata

export seriesoverlay
"""
    seriesoverlay(ts1, ts2)

Return a new [`TSeries`](@ref) over the full range of both arguments. The overlapping part
contains values from the last argument.

See also: [`dictoverlay`](@ref)

"""
function seriesoverlay(ts1::TSeries, ts2::TSeries)
    # Make a copy of the output sries
    tsout = deepcopy(ts2);
    # Range of the first TSeries
    Rng1 = mitrange(ts1);
    # Range of the second TSeries
    Rng2 = mitrange(ts2);
    if Rng1.start < Rng2.start
        tsout[Rng1.start:Rng2.start - 1] = ts1;
    end
    if Rng1.stop > Rng2.stop
        tsout[Rng2.stop + 1:Rng1.stop] = ts1;
    end
    return tsout
end

export dictoverlay
"""
    dictoverlay(d1, d2)

Merge two dictionaries. Common key where the values are [`TSeries`](@ref) of the
same frequency are overlayed. Otherwise, a common key takes the value of the
last Dict containing it.

See also: [`seriesoverlay`](@ref)
"""
function dictoverlay(D1::Dict{String,<:Any}, D2::Dict{String,<:Any})
    # Get keys
    K1 = keys(D1);
    K2 = keys(D2);
    # Combine the keys
    K = unique([collect(K1) ; collect(K2)]);
    # Pre-allocate output
    D3 = Dict{String,Any}();
    # Go over each key and update D3
    for k in K
        # Test where we can find the key
        t1 = in(k, K1)
        t2 = in(k, K2)
        if !t1
            # If we cannot find it in K1, we use D2
            push!(D3, k => D2[k])
        elseif !t2
            # If we cannot find it in K2, we use D1
            push!(D3, k => D1[k])
        else
            # If we find it in both D1 and D2
            if isa(D1[k], TSeries) && isa(D2[k], TSeries)
                # If both are TSeries, we overlay them.
                push!(D3, k => seriesoverlay(D1[k], D2[k]))
            else
                # Otherwise, we give priority to D2
                push!(D3, k => D2[k])
            end
        end
    end
    return D3 = deepcopy(D3)
end

export array2dict
"""
    array2dict(data, vars, start_date)

Convert the simulation data array to a dictionary.
"""
function array2dict(data::AbstractArray{Float64,2}, vars::AbstractVector, start_date::MIT)::Dict{String,Any}
    Dict{String,Any}(string(vars[i]) => TSeries(start_date, data[:,i]) for i ∈ 1:length(vars))
end

export array2data
"""
    array2data(data, vars, start_date; copy=false)

Convert the simulation data array to a [`SimData`](@ref).

Use the `copy` argument to control whether or not the returned SimData will hold
a reference to the given data array or its own copy.

"""
function array2data(data::AbstractArray{Float64,2}, vars::AbstractVector, start_date::MIT; copy=false)::SimData
    SimData(start_date, vars, ifelse(copy, Base.copy(data), data))
end

export dict2array

function _d2_vars(d, vars)
    if keytype(d) <: AbstractString
        vars = string.(vars)
    else
        vars = keytype(d).(vars)
    end
    missing_vars = setdiff(vars, keys(d))
    if !isempty(missing_vars)
        throw(ArgumentError("""Variables not found in data dictionary: $(join(missing_vars, ", "))"""))
    end
    return vars
end

function _d2_rng(d, vars, range)
    ranges = [mitrange(d[v]) for v in vars]
    fs = first.(ranges)
    ls = last.(ranges)
    if isempty(range)
        f = maximum(fs)
        l = minimum(ls)
        range = f:l
        if !(all(f .== fs) && all(l .== ls))
            @warn "Variable ranges are not all the same. Specify `range=` to avoid this warning. Using intersection range: $(range)"
        end
    else
        missing_data = (first(range) .< fs) .| (last(range) .> ls)
        if any(missing_data)
            throw(ArgumentError("""Data is not available for the full range $(range) for some variables: $(join(vars[missing_data], ", "))"""))
        end
    end
    return range
end

"""
    dict2array(d, vars; range)

Convert a dictionary of [`TSeries`](@ref) to a 2d array of simulation data for
the given range. If the `range` argument is not provided, the effective range is
the intersection of the ranges of available data for the given list of
variables.

"""
function dict2array(d::Dict, vars::AbstractVector; range::AbstractUnitRange=1:0)::Array{Float64,2}
    vars = _d2_vars(d, vars)
    range = _d2_rng(d, vars, range)
    data = Array{Float64,2}(undef, length(range), length(vars))
    for i in eachindex(vars)
        data[:, i] .= d[vars[i]][range].values
    end
    return data
end

export dict2data
"""
    dict2data(d, vars; range)

Convert a dictionary of [`TSeries`](@ref) to a `SimData`` for the given range.
If the `range` argument is not provided, the effective range is the intersection
of the ranges of available data for the given list of variables.

"""
function dict2data(d::Dict, vars::AbstractVector; range::AbstractUnitRange=1:0)::SimData
    vars = _d2_vars(d, vars)
    range = _d2_rng(d, vars, range)
    data = Array{Float64,2}(undef, length(range), length(vars))
    for i in eachindex(vars)
        data[:, i] .= d[vars[i]][range].values
    end
    return SimData(first(range), vars, data)
end

export data2dict
"""
    data2dict(sd::SimData; copy=false)::Dict{String, Any}

Convert a [`SimData`](@ref) to a dictionary containing the same data as
individual [`TSeries`](@ref).

Use the `copy` argument to control whether the TSeries in the returned Dict will
hold references to the columns of `sd` or copies of that data.

"""
function data2dict(sd::SimData; copy=false)
    transform = ifelse(copy, Base.copy, Base.identity)
    Dict{String,Any}(string(n) => transform(v) for (n, v) in pairs(sd))
end

export data2array
"""
    data2array(sd::SimData; copy=false)::Array{Float64,2}

Convert a [`SimData`](@ref) to an Array. The `copy` argument controls whether
the returned Array holds a copy of the data or a reference to the data in `sd`.

"""
function data2array(sd::SimData; copy=false)
    return ifelse(copy, Base.copy(rawdata(sd)), rawdata(sd))
end


text:
mutable struct Obstacle
    body::RigidBody
    point::Point3D
    normal::FreeVector3D
    basis::Matrix{Float64}
    μ::Real
    is_floating::Bool

    function Obstacle(body::RigidBody, point::Point3D, normal::FreeVector3D, motion_type::Symbol, μ::Real; is_floating=false)
        @framecheck point.frame normal.frame
        normal_norm = normalize(normal)
        basis = contact_basis(normal_norm, μ, motion_type)
        new(body, point, normal_norm, basis, μ, is_floating)
    end
end

separation(obs::Obstacle, p::Point3D) = dot(p - obs.point, obs.normal)

ext: jl


text:
function test_detector(detector::OD.Detector, data::TestData)
    @testset "$detector" begin
        is_supervised = isa(detector, OD.SupervisedDetector)

        # OutlierDetection detector with matrix input
        raw_model, raw_train = is_supervised ?
            OD.fit(detector, data.x_raw[:, data.train_idx], data.y_raw[data.train_idx]; verbosity=0) :
            OD.fit(detector, data.x_raw[:, data.train_idx]; verbosity=0)
        raw_test = OD.transform(detector, raw_model, data.x_raw[:, data.test_idx])

        # MLJ plain detector with table input
        mlj_raw = is_supervised ? MLJ.machine(detector, data.x_table, data.y_cat) :
                                  MLJ.machine(detector, data.x_table)
        MLJ.fit!(mlj_raw, rows=data.train_idx)
        mlj_raw_train = MLJ.report(mlj_raw).scores
        mlj_raw_test = MLJ.transform(mlj_raw, rows=data.test_idx)

        # MLJ probabilistic detector with table input
        mlj_prob = is_supervised ? MLJ.machine(ProbabilisticDetector(detector), data.x_table, data.y_cat) :
                                   MLJ.machine(ProbabilisticDetector(detector), data.x_table)
        MLJ.fit!(mlj_prob, rows=data.train_idx)
        mlj_prob_train = MLJ.report(mlj_prob).detector.scores
        mlj_prob_test = from_univariate_finite(MLJ.predict(mlj_prob, rows=data.test_idx))

        # MLJ deterministic detector with table input
        mlj_det = is_supervised ? MLJ.machine(DeterministicDetector(detector), data.x_table, data.y_cat) :
                                  MLJ.machine(DeterministicDetector(detector), data.x_table)
        MLJ.fit!(mlj_det, rows=data.train_idx)
        mlj_det_train = MLJ.report(mlj_det).detector.scores
        mlj_det_test = from_categorical(MLJ.predict(mlj_det, rows=data.test_idx))

        @testset "outputs have appropriate dimensions" begin
            @test length(raw_train) ==
                  length(mlj_raw_train) ==
                  length(mlj_prob_train) ==
                  length(mlj_det_train) ==
                  data.n_train

            @test length(raw_test) ==
                  length(mlj_raw_test) ==
                  length(mlj_prob_test) ==
                  length(mlj_det_test) ==
                  data.n_test
        end

        @testset "raw scores have appropriate values" begin
            @test all(-Inf .< raw_train .< Inf)
            @test all(-Inf .< mlj_raw_train .< Inf)
            @test all(-Inf .< mlj_prob_train .< Inf)
            @test all(-Inf .< mlj_det_train .< Inf)
            @test all(-Inf .< raw_test .< Inf)
            @test all(-Inf .< mlj_raw_test .< Inf)
        end

        @testset "normalized scores have appropriate values" begin
            @test all(0 .<= mlj_prob_test .<= 1)
        end

        @testset "labels have appropriate values" begin
            in_labels(x) = x in (OutlierDetection.CLASS_NORMAL, OutlierDetection.CLASS_OUTLIER)
            @test all(in_labels.(mlj_det_test))
        end
    end
end


text:
export circularShutter!, circularShutterFreq!, getCircularMask

preserveType(orig::Array{T}, data::Array{U}) where {T<:Real, U<:Complex} = real(data)

preserveType(orig::Array{T}, data::Array{U}) where {T, U} = data

function circularShutter!(I::Matrix, radiusFactor::Number=1.0)
    imgSize = size(I)
    center = (imgSize[1]/2.0, imgSize[2]/2.0)
    radius = (center[1]>center[2] ? center[1] : center[2]) * radiusFactor
    # Applying filtering
    for j=1:imgSize[2]
        for i=1:imgSize[1]
            if sqrt((i-center[1])^2 + (j-center[2])^2) > radius
                I[i,j] = 0
            end
        end
    end
end

function circularShutter!(I::Array{Complex{T}, 3}, radiusFactor::Float64=1.0) where T <: AbstractFloat
    imgSize = size(I)
    center = (imgSize[1]/2.0, imgSize[2]/2.0, imgSize[3]/2.0)
    radius = maximum(center) * radiusFactor
    # Applying filtering
    for k=1:imgSize[3]
        for j=1:imgSize[2]
            for i=1:imgSize[1]
                if sqrt((i-center[1])^2 + (j-center[2])^2 + (k-center[3])^2) > radius
                    I[i,j,k] = 0
                end
            end
        end
    end
end

function getCircularMask(shape::Tuple, radiusFactor::Number=1.0)
  A = ones(Float64,shape)
  circularShutter!(A,radiusFactor)
  return A
end

function circularShutterFreq!(I::Matrix, radiusFactor::T=1.0) where T<:Number
    imgSize = size(I)
    fftI = fftshift(fft(I))
    center = (imgSize[1]/2.0, imgSize[2]/2.0)
    radius = (center[1]>center[2] ? center[1] : center[2]) * radiusFactor
    # Applying filtering
    for j=1:imgSize[2]
        for i=1:imgSize[1]
            if sqrt((i-center[1])^2 + (j-center[2])^2) > radius
                fftI[i,j] = 0
            end
        end
    end
    return preserveType(I, ifft(ifftshift(fftI)))
end

function circularShutterFreq!(I::Array{T,3}, radiusFactor::T=1.0) where T<:Number
    imgSize = size(I)
    fftI = fftshift(fft(I))
    center = (imgSize[1]/2.0, imgSize[2]/2.0, imgSize[3]/2.0)
    radius = (center[1]>center[2] ? center[1] : center[2]) * radiusFactor
    # Applying filtering
    for k=1:imgSize[3]
        for j=1:imgSize[2]
          for i=1:imgSize[1]
            if sqrt((i-center[1])^2 + (j-center[2])^2 + (k-center[3])^2) > radius
                fftI[i,j,k] = 0
            end
          end
        end
    end
    return preserveType(I, ifft(ifftshift(fftI)))
end


text:
# Tuple Tests
let
    basetup = (1,2.0,3,4.0,5,6.0)
    for j in 1:length(basetup)
        tup = tuple(basetup[1:j])
        for i in 1:length(tup)
            @test typemax(tup)[i] == typemax(typeof(tup[i]))
            @test typemin(tup)[i] == typemin(typeof(tup[i]))
        end
    end
    @test (1,0,0,0,0) < 2
    @test (1,0,0,0,0) <= 1
    @test (1,0,0,0,0) >= 1
    @test (1,0,0,0,0) > 0

    @test 2 > (1,0,0,0,0)
    @test 1 >= (1,0,0,0,0)
    @test 1 <= (1,0,0,0,0)
    @test 0 < (1,0,0,0,0)
end
let
    S = GraphState
    A = GraphAction
    P = PathNode{S,A}

    Path([P()])
    p = Path{S,A,Float64}()
    @test state_type(p) == S
    @test action_type(p) == A
    @test node_type(p) == P
    @test cost_type(p) == Float64
    get(p,1,node_type(p))
    get(p,1)
    # @test get_cost(p) == 0
    @test length(p.path_nodes) == 0
    push!(p,P())
    @test length(p.path_nodes) == 1
    p = cat(p,P())
    @test length(p.path_nodes) == 2
    @test typeof(p[1]) == P
    p[1] = P()
    @test length(p) == length(p.path_nodes)
    p1 = copy(p)
    @test get_cost(p1) == get_cost(p)
    @test length(p1) == length(p)
end
let
    S = GraphState
    A = GraphAction
    P = PathNode{S,A}

    p = Path{S,A,Float64}(
        s0 = S(1,0)
    )

    get_path_node(p,1)
    get_path_node(p,5)
    get_a(p,1)
    @test length(extend_path(p,5)) == 5
    @test get_end_index(extend_path(p,5)) == 5
    extend_path!(p,5)
    @test length(p) == 5
    @test get_end_index(p) == 5

    get_initial_state(p)
    get_final_state(p)

    cost_model = TravelTime()
    solution = LowLevelSolution(
        cost_model = TravelTime(),
        paths = [Path{S,A,cost_type(cost_model)}()],
        costs = [0.0]
    )
    state_type(solution)
    action_type(solution)
    cost_type(solution)
    node_type(solution)
    copy(solution)
    get_paths(solution)
    get_path_costs(solution)
    get_cost(solution)
    get_cost_model(solution)
    set_solution_path!(solution,get_paths(solution)[1],1)
    set_path_cost!(solution,2.0,1)
end


text:
## Framework for setting up an iterative problem for finding a zero
## Different methods can implement:
## - A subtype of UnivariateZeroMethod (reqd)
## - UnivariateZeroProblem
## - callable_function
## - init_state
## - assess_convergence
## - udpate_state (reqd)
## TODO
## * a graphic of trace when verbose=true?

## method names are subtypes
abstract type UnivariateZeroMethod end
abstract type AbstractBisection <: UnivariateZeroMethod end
abstract type AbstractSecant <: UnivariateZeroMethod end

# container for callable objects; not really necessary, but has some value.
abstract type CallableFunction end
struct DerivativeFree{F} <: CallableFunction
    f::F
end
(F::DerivativeFree)(x::Number) = F.f(x)

struct FirstDerivative{F,D} <: CallableFunction
    f::F
    fp::D
end
(F::FirstDerivative)(x::Number) = F.f(x)

struct SecondDerivative{F,D,H} <: CallableFunction
    f::F
    fp::D
    fpp::H
end
(F::SecondDerivative)(x::Number) = F.f(x)

## allows override for automatic derivatives, see Newton
function callable_function(m::UnivariateZeroMethod, @nospecialize(f))
    !isa(f, Tuple) && return DerivativeFree(f)
    length(f) == 1 && return DerivativeFree(f[1])
    length(f) == 2 && return FirstDerivative(f[1], f[2])
    SecondDerivative(f[1], f[2], f[3])
end

## object to hold state; allow for subtypes with additional fields
abstract type  UnivariateZeroState{T, S} end
mutable struct UnivariateZeroStateBase{T,S} <: UnivariateZeroState{T,S}
    xn1::T
    xn0::T
    fxn1::S
    fxn0::S
    bracket::Union{Vector{T}, Missing}
    steps::Int
    fnevals::Int
    stopped::Bool             # stopped, butmay not have converged
    x_converged::Bool         # converged via |x_n - x_{n-1}| < ϵ
    f_converged::Bool         # converged via |f(x_n)| < ϵ
    convergence_failed::Bool
    message::AbstractString
end

incfn(o::UnivariateZeroState, k=1)    = o.fnevals += k
incsteps(o::UnivariateZeroState, k=1) = o.steps += k

## generic initialization. Modify as necessary for a method, such as secant which uses both xn1 and xn0
## we use x0 + typemax(Int) as a sentinel. This could be a Missing, but that
## is a bit more hassle
function init_state(method::Any, fs, x0::T, bracket) where {T}
    fx0 = fs(x0); fnevals = 1

    if !ismissing(bracket)
        a,b = bracket[1:2]
        sign(fs(a)) * sign(fs(b)) > 0 &&  (warn(bracketing_error); throw(ArgumentError))
        fnevals += 2
    end

    S = eltype(fx0)
    state = UnivariateZeroStateBase(x0, x0 + typemax(Int),
                                    fx0, fx0,
                                    !ismissing(bracket) ? float.(bracket) : bracket,
                                    0, fnevals,
                                    false, false, false, false,
                                    "")
    state
end

## Options for convergence, reporting
mutable struct UnivariateZeroOptions{T}
    xabstol::T
    xreltol::T
    abstol::T
    reltol::T
    maxevals::Int
    maxfnevals::Int
    verbose::Bool
end

function univariate_zero_options(args...;
                                    xabstol::T=zero(T),
                                    xreltol::T=eps(T),
                                    abstol::T=zero(T),
                                    reltol::T=4*eps(T),
                                    maxevals::Int=40,
                                    maxfnevals::Int=typemax(Int),
                                    verbose::Bool=false,
                                    kwargs...) where {T}

    ## adjust for old argument names
    kw = Dict(kwargs)
    z = zero(T)
    UnivariateZeroOptions(max(z, get(kw, :xtol, xabstol)),
                          max(z, get(kw, :xtolrel, xreltol)),
                          max(z, abstol),
                          max(z, get(kw, :ftol, reltol)),
                          max(0, get(kw, :maxeval, get(kw, :maxsteps, maxevals))),
                          max(0, get(kw, :maxfneval, maxfnevals)),
                          verbose)
end

## basic container
mutable struct UnivariateZeroProblem{T<:AbstractFloat, F, XT <: Union{T, Tuple{T,T}, Vector{T}, Complex{T}}, B <: Union{Vector{T}, Missing}}
    fs::F
    x0::XT
    bracket::B
end
UnivariateZeroProblem(fs::F, x0::XT, bracket::Vector{T}) where {T,F,XT<:Union{T,Tuple{T,T},Vector{T},Complex{T}}} = UnivariateZeroProblem{T,F,XT,Vector{T}}(fs,x0,bracket)
UnivariateZeroProblem(fs::F, x0::T, bracket::Missing) where {T<:AbstractFloat,F} = UnivariateZeroProblem{T,F,T,Missing}(fs,x0,bracket)
UnivariateZeroProblem(fs::F, x0::Tuple{T,T}, bracket::Missing) where {T<:AbstractFloat,F} = UnivariateZeroProblem{T,F,Tuple{T,T},Missing}(fs,x0,bracket)
UnivariateZeroProblem(fs::F, x0::Vector{T}, bracket::Missing) where {T<:AbstractFloat,F} = UnivariateZeroProblem{T,F,Vector{T},Missing}(fs,x0,bracket)
UnivariateZeroProblem(fs::F, x0::Complex{T}, bracket::Missing) where {T<:AbstractFloat,F} = UnivariateZeroProblem{T,F,Complex{T},Missing}(fs,x0,bracket)

## frame the problem and the options
function derivative_free_setup(method::Any, fs, x0::Union{T, Tuple{T,T}, Vector{T}};
                               bracket=missing,
                               xabstol=zero(T), xreltol=eps(T),
                               abstol=4*eps(T), reltol=4*eps(T),
                               maxevals=40, maxfnevals=typemax(Int),
                               verbose::Bool=false,
                               kwargs...
    ) where {T<:AbstractFloat}
    x = float.(x0)
    prob = UnivariateZeroProblem(fs, x,
                                 !ismissing(bracket) ? float.(bracket) : missing)
    options = UnivariateZeroOptions(xabstol, xreltol, abstol,
                                    reltol, maxevals, maxfnevals, verbose)
    prob, options
end

## has UnivariateZeroProblem converged?
function assess_convergence(method::Any, fs, state, options)

    xn0, xn1 = state.xn0, state.xn1
    fxn0, fxn1 = state.fxn0, state.fxn1

    if (state.x_converged || state.f_converged)
        return true
    end

    if state.steps > options.maxevals
        state.stopped = true
        state.message = "too many steps taken."
        return true
    end

    if state.fnevals > options.maxfnevals
        state.stopped = true
        state.message = "too many function evaluations taken."
        return true
    end

    if isnan(xn1)
        state.convergence_failed = true
        state.message = "NaN produced by algorithm"
        return true
    end

    if isinf(fxn1)
        state.convergence_failed = true
        state.message = "Inf produced by algorithm"
        return true
    end

    λ = max(one(real(xn1)), norm(xn1))

    if  norm(fxn1) <= max(options.abstol, λ * options.reltol)
        state.f_converged = true
        return true
    end

    if isapprox(xn1, xn0, rtol = options.xreltol, atol=options.xabstol) && norm(fxn1) <= cbrt(max(options.abstol, λ * options.reltol))
        state.x_converged = true
        return true
    end

    if state.stopped
        if state.message == ""
            error("no message? XXX debug this XXX")
        end
        return true
    end

    return false
end

function show_trace(fs, state, xns, fxns, method)
    converged = state.x_converged || state.f_converged

    println("Results of univariate zero finding:\n")
    if converged
        println("* Converged to: $(xns[end])")
        println("* Algorithm: $(method)")
        println("* iterations: $(state.steps)")
        println("* function evaluations: $(state.fnevals)")
        state.x_converged && println("* stopped as x_n ≈ x_{n-1} using atol=xabstol, rtol=xreltol")
        state.f_converged && state.message == "" && println("* stopped as |f(x_n)| ≤ max(δ, max(1,|x|)⋅ϵ) using δ = abstol, ϵ = reltol")
        state.message != "" && println("* Note: $(state.message)")
    else
        println("* Convergence failed: $(state.message)")
        println("* Algorithm $(method)")
    end
    println("")
    println("Trace:")

    itr, offset =  0:(endof(xns)-1), 1
    for i in itr
        x_i,fx_i, xi, fxi = "x_$i", "f(x_$i)", xns[i+offset], fxns[i+offset]
        println(@sprintf("%s = % 18.16f,\t %s = % 18.16f", x_i, float(xi), fx_i, float(fxi)))
    end
    println("")

end

### find_zero method has potentially many interfaces

## problem, method, options. This may be a main entry point
function find_zero(prob::UnivariateZeroProblem,  method::UnivariateZeroMethod, options::UnivariateZeroOptions)
    fs = prob.fs
    state = init_state(method, fs, prob.x0, prob.bracket)
    find_zero(method, fs, state, options)
end

## method, state, options. Could be used to switch  between methods
function find_zero(method::UnivariateZeroMethod, fs, state::UnivariateZeroState, options::UnivariateZeroOptions)

    # in case verbose=true
    if isa(method, AbstractSecant)
        xns, fxns = [state.xn0, state.xn1], [state.fxn0, state.fxn1]
    else
        xns, fxns = [state.xn1], [state.fxn1]
    end

    ## XXX Should just deprecate this in favor of FalsePosition method XXX
    while true
        if  !ismissing(state.bracket)
            m,M = state.bracket[1:2]
            if (state.xn1 < m || state.xn1 > M) || state.stopped
                # do bisection step, update bracket
                c = m + (M-m)/2
                # don't land on xn0
                if c == state.xn0
                    c = m + (M-m)/4
                end
                fc = fs(c)
                if norm(fc) < norm(state.fxn1)
                    state.xn0 = state.xn1
                    state.xn1 = c
                    state.fxn0 = state.fxn1
                    state.fxn1 = fc
                else
                    state.xn0 = c
                    state.fxn1 = fc
                end

                if sign(state.fxn1) * sign(fs(m)) < 0
                    state.bracket[2] = c
                else
                    state.bracket[1] = c
                end
                state.stopped && (state.stopped = false)
                incfn(state)
            end
        end

        val = assess_convergence(method, fs, state, options)

        if val
            if state.stopped && !(state.x_converged || state.f_converged)
                ## stopped is a heuristic, there was an issue with an approximate derivative
                ## say it converged if pretty close, else say convergence failed.
                ## (Is this a good idea?)
                xstar, fxstar = state.xn1, state.fxn1
                if abs(fxstar) <= (options.abstol)^(2/3)
                    msg = "Algorithm stopped early, but |f(xn)| < ϵ^(2/3), where ϵ = abstol"
                    state.message = state.message == "" ? msg : state.message * "\n\t" * msg
                    state.f_converged = true
                else
                    state.convergence_failed = true
                end
            end

            if state.x_converged || state.f_converged
                options.verbose && show_trace(fs, state, xns, fxns, method)
                return state.xn1
            end

            if state.convergence_failed
                options.verbose && show_trace(fs, state, xns, fxns, method)
                throw(ConvergenceFailed("Stopped at: xn = $(state.xn1)"))
            end
        end

        update_state(method, fs, state, options)
        if options.verbose
            push!(xns, state.xn1)
            push!(fxns, state.fxn1)
        end

    end
end

"""

Find a zero of a univariate function using one of several different methods.

Positional arugments:

* `f` a function, callable object, or tuple of same. A tuple is used
  to pass in derivatives, as desired. Most methods are derivative
  free. Some (`Newton`, `Halley`) may have derivative(s) computed
  using the `ForwardDiff` pacakge.

* `x0` an initial starting value. Typically a scalar, but may be a two-element
  tuple or array for bisection methods. The value `float.(x0)` is passed on.

* `method` one of several methods, see below.

Keyword arguments:

* `xabstol=zero()`: declare convergence if |x_n - x_{n-1}| <= max(xabstol, max(1, |x_n|) * xreltol)

* `xreltol=eps()`:

* `abstol=zero()`: declare convergence if |f(x_n)| <= max(abstol, max(1, |x_n|) * reltol)

* `reltol`:

* `bracket`: Optional. A bracketing interval for the sought after
  root. If given, a hybrid algorithm may be used where bisection is
  utilized for steps that would go out of bounds. (Using a `FalsePosition` method
  instead would be suggested.)

* `maxevals::Int=40`: stop trying after `maxevals` steps

* `maxfnevals::Int=typemax(Int)`: stop trying after `maxfnevals` function evaluations

* `verbose::Bool=false`: If `true` show information about algorithm and a trace.

Returns:

Returns `xn` if the algorithm converges. If the algorithm stops, returns `xn` if
|f(xn)| ≤ ϵ^(2/3), where ϵ = reltol, otherwise a `ConvergenceFailed` error is thrown.

Exported methods:

`Bisection()`;
`Order0()` (heuristic, slow more robust);
`Order1()` (also `Secant()`);
`Order2()` (also `Steffensen()`);
`Order5()` (KSS);
`Order8()` (Thukral);
`Order16()` (Thukral);
`FalsePosition(i)` (false position, i in 1..12);

Not exported:

`Secant()`, use `Order1()`
`Steffensen()` use `Order2()`
`Newton()` (use `newton()` function)
`Halley()` (use `halley()` function)

The order 0 method is more robust to the initial starting point, but
can utilize many more function calls. The higher order methods may be
of use when greater precision is desired.`

Examples:

```
f(x) = x^5 - x - 1
find_zero(f, 1.0, Order5())
find_zero(f, 1.0, Steffensen()) # also Order2()
find_zero(f, (1.0, 2.0), FalsePosition())
```
"""
function find_zero(f, x0::T, method::M; kwargs...) where {M <: AbstractBisection, T<:Number}
    throw(ArgumentError("For bisection methods, x0 must be a vector"))
end

function find_zero(f, x0::Vector{T}, method::M; kwargs...) where {T<:Number, M<:AbstractBisection}
    find_zero(method, callable_function(method, f), x0; kwargs...)
end

function find_zero(f, x0::Tuple{T}, method::M; kwargs...) where {T<:Number, M<:AbstractBisection}
    find_zero(method, callable_function(method, f), x0; kwargs...)
end

function find_zero(f, x0::T, method::UnivariateZeroMethod; kwargs...) where {T<:Number}
    find_zero(method, callable_function(method, f), x0; kwargs...)
end

function find_zero(f, x0::T, method::AbstractSecant; kwargs...) where {T<:Number}
    find_zero(method, callable_function(method, f), x0; kwargs...)
end
function find_zero(f, x0::Vector{T}, method::AbstractSecant; kwargs...) where {T<:Number}
    find_zero(method, callable_function(method, f), x0; kwargs...)
end

## some defaults for methods
find_zero(f, x0::T; kwargs...)  where {T <: Number}= find_zero(f, x0, Order0(); kwargs...)
find_zero(f, x0::Vector{T}; kwargs...) where {T <: Number}= find_zero(f, x0, Bisection(); kwargs...)
find_zero(f, x0::Tuple{T,S};kwargs...) where {T<:Number, S<:Number} = find_zero(f, x0, Bisection(); kwargs...)

function find_zero(method::UnivariateZeroMethod, fs, x0; kwargs...)
    x = float.(x0)
    prob, options = derivative_free_setup(method, fs, x; kwargs...)
    find_zero(prob, method, options)
end

## old interface for fzero
## old keyword arguments (see ?fzero) handled in univariate_zero_options
@noinline function derivative_free(f, x0::T; order::Int=0,
                                             kwargs...) where {T <: AbstractFloat}

    if order == 0
        method = Order0()
    elseif order == 1
        method = Order1()
    elseif order == 2
        method = Order2()
    elseif order == 5
        method = Order5()
    elseif order == 8
        method = Order8()
    elseif order == 16
        method = Order16()
    else
        warn("Invalid order. Valid orders are 0, 1, 2, 5, 8, and 16")
        throw(ArgumentError())
    end

    find_zero(f, x0, method; kwargs...)
end


text:
## Some more extreme tests from http://ir.igsnrr.ac.cn/bitstream/311030/8840/1/%E4%BE%AF%E9%BA%9F%E7%A7%91(SCI)2.pdf
## All the methods do poorly here due to the multiplicities of the functions
## even though there should be convergnce of either |xn+1 - xn| or |f(xn)|
using Roots
using Base.Test

tests = [
         (x -> (x - sqrt(5))^4 / ((x-1)^2 + 2),    4.3,  2.236067977499790)
         (x -> (sin(x)^2 - 2x + 1)^5,              4.5,  0.71483582544138924)
         (x -> (8x*exp(-x^2) -2x - 3)^8,          -2.0, -1.7903531791589544)
         (x -> (2x*cos(x) + x^2 - 3)^10/(x^2 + 1), 3.0,  2.9806452794385368)
         (x -> (exp(-x^2 + x + 3) - x + 2)^9,      2.5,  2.4905398276083051)
         (x -> (exp(-x) + 2sin(x))^4,              3.0,  3.1627488709263654)
         (x -> (log(x^2 + 3x + 5) - 2x + 7)^8,     5.5,  5.4690123359101421)
         (x -> (sqrt(x^2 + 2x + 5) -
                           2sin(x) - x^2 + 3)^5,   2.3, 2.3319676558839640)
         (x -> (x-2)^4 / ( (x-1)^2 + 1),           1.8,  2.0000000000000000)
         (x -> abs(x - 2.5)^(15/4) * exp(x),       2.2,  2.5)
         (x -> (sqrt(x) - 1/x - 1)^7,              2.0,  2.147899035704787)
         (x -> (log(x) + sqrt(x) - 5)^3,           8.0,  8.309432694231572)
         (x -> (sin(x)*cos(x) - x^3 + 1)^9,        1.0,  1.117078770687451)
         (x -> ((x-3)*exp(x))^5,                   2.5,  3.0000000000000000)
         (x -> (log(x) + sqrt(x^4 + 1) -2)^7,      1.0,  1.222813963628973)
         ]

using Roots
ord(x) = string(floor(Integer,log10(x)))
orders = [0, 1, 2, 5, 8, 16]
out = Any[]
for (ctr, (f1, x0, xstar)) in enumerate(tests)
    m = AbstractString[]
    for order in orders
        try
            a = fzero(f1, x0, order=order, maxeval=100)
            tm = begin tic(); fzero(f1, x0, order=order, maxeval=100); toc() end
            push!(m, ord(tm))
            push!(m, ord(abs(a-xstar)))
        catch
            push!(m, "***")
            push!(m, "***")
        end
    end
    try
        newton(f1, x0)
        tm =  begin tic(); a = newton(f1, x0); toc() end
        push!(m, ord(tm))
        push!(m, ord(abs(a-xstar)))
    catch e
        push!(m, "***")
        push!(m, "***")
    end
    push!(out, m)
end

## values in pairs time/accuracy or ****** if failure
## sensitive to maxeval value
 for i in out
     for k in 1:div(length(i),2)
         @printf "%3s%3s " i[2*(k-1)+1] i[2*k]
     end
         @printf "\n"
 end


text:
include("plot_match.jl")

@recipe function plot(avg_wave::AverageWave{T};
        hankel_indexes = 0:avg_wave.hankel_order,
        apply = real) where T<:AbstractFloat

    ho = avg_wave.hankel_order

    for n in hankel_indexes

        apply_field = apply.(avg_wave.amplitudes[:,n+ho+1,1])

        @series begin
            label --> "$apply Hankel = $n"
            seriestype --> :scatter
            (avg_wave.x, apply_field)
        end
    end
end

@recipe function plot(wave_effs::Vector{EffectiveWave{T}}) where T<:AbstractFloat
    k_effs = [w.k_eff for w in wave_effs]
    maxamp = maximum(norm(w.amplitudes) for w in wave_effs)

    alphas = map(wave_effs) do w
        norm(w.amplitudes)/maxamp
    end

    @series begin
        ylims --> (0,Inf)
        xlab --> "Re k_eff"
        ylab --> "Im k_eff"
        seriestype --> :scatter
        label --> ""
        markerstrokealpha --> 0.2
        # markercolor --> :blue
        markeralpha --> alphas
        k_effs
    end
end

@recipe function plot(x::AbstractVector{T}, wave_eff::EffectiveWave{T}) where T<:AbstractFloat
    @series begin
        (x, [wave_eff])
    end
end

@recipe function plot(x::AbstractVector, wave_effs::Vector{E};
        hankel_indexes = 0:wave_effs[1].hankel_order,
        apply = real) where E<:EffectiveWave

    wave_eff = AverageWave(x, wave_effs)
    ho = wave_eff.hankel_order

    for n in hankel_indexes

        apply_field = apply.(wave_eff.amplitudes[:,n+ho+1,1])

        @series begin
            label --> "$apply Hankel = $n"
            (x, apply_field)
        end
    end
end


text:
#! format: off

const ACTIVE_POWER = "P"
const ACTIVE_POWER_IN = "Pin"
const ACTIVE_POWER_OUT = "Pout"
const COLD_START = "start_cold"
const ENERGY = "E"
const ENERGY_UP = "Eup"
const ENERGY_DOWN = "Edown"
const ENERGY_BUDGET = "energy_budget"
const ENERGY_BUDGET_UP = "energy_budget_up"
const ENERGY_BUDGET_DOWN = "energy_budget_down"
const ENERGY_SHORTAGE  = "energy_shortage"
const ENERGY_SURPLUS = "energy_surplus"
const FLOW_REACTIVE_POWER_FROM_TO = "FqFT"
const FLOW_REACTIVE_POWER_TO_FROM = "FqTF"
const FLOW_ACTIVE_POWER_FROM_TO = "FpFT"
const FLOW_ACTIVE_POWER_TO_FROM = "FpTF"
const FLOW_ACTIVE_POWER = "Fp"
const FLOW_REACTIVE_POWER = "Fq"
const HOT_START = "start_hot"
const INFLOW = "In"
const TARGET = "Target"
const OUTFLOW = "Out"
const ON = "On"
const REACTIVE_POWER = "Q"
const RESERVE = "R"
const SERVICE_REQUIREMENT = "service_requirement"
const SLACK_DN = "γ⁻"
const SLACK_UP = "γ⁺"
const SPILLAGE = "Sp"
const START = "start"
const STOP = "stop"
const THETA = "theta"
const VM = "Vm"
const LIFT = "z"
const ACTIVE_POWER_PUMP = "Ppump"

abstract type VariableType end

"""Struct to dispatch the creation of Active Power Variables"""
struct ActivePowerVariable <: VariableType end

"""Struct to dispatch the creation of Active Power Input Variables for 2-directional devices. for instance storage or pump-hydro"""
struct ActivePowerInVariable <: VariableType end

"""Struct to dispatch the creation of Active Power Output Variables for 2-directional devices. for instance storage or pump-hydro"""
struct ActivePowerOutVariable <: VariableType end

struct HotStartVariable <: VariableType end

struct WarmStartVariable <: VariableType end

struct ColdStartVariable <: VariableType end

struct EnergyVariable <: VariableType end

struct EnergyVariableUp <: VariableType end

struct EnergyVariableDown <: VariableType end

struct EnergyShortageVariable <: VariableType end

struct EnergySurplusVariable <: VariableType end

struct LiftVariable <: VariableType end

struct OnVariable <: VariableType end

struct ReactivePowerVariable <: VariableType end

struct ReserveVariable <: VariableType end

struct ActiveServiceVariable <: VariableType end

struct ServiceRequirementVariable <: VariableType end

struct SpillageVariable <: VariableType end

struct StartVariable <: VariableType end

struct StopVariable <: VariableType end

struct SteadyStateFrequencyDeviation <: VariableType end

struct AreaMismatchVariable <: VariableType end

struct DeltaActivePowerUpVariable <: VariableType end

struct DeltaActivePowerDownVariable <: VariableType end

struct AdditionalDeltaActivePowerUpVariable <: VariableType end

struct AdditionalDeltaActivePowerDownVariable <: VariableType end

struct SmoothACE <: VariableType end

"""Struct to dispatch the creation of Flow Active Power Variables"""
struct FlowActivePowerVariable <: VariableType end

struct FlowReactivePowerVariable <: VariableType end

struct FlowActivePowerFromToVariable <: VariableType end

struct FlowActivePowerToFromVariable <: VariableType end

struct FlowReactivePowerFromToVariable <: VariableType end

struct FlowReactivePowerToFromVariable <: VariableType end

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

const START_VARIABLES = (HotStartVariable, WarmStartVariable, ColdStartVariable)

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

make_variable_name(var_type, device_type) = encode_symbol(device_type, var_type)
make_variable_name(var_type) = encode_symbol(var_type)

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

make_variable_name(::Type{ActivePowerVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "P")

make_variable_name(::Type{ActivePowerInVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "Pin")

make_variable_name(::Type{ActivePowerOutVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "Pout")

make_variable_name(::Type{HotStartVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "start_hot")

make_variable_name(::Type{WarmStartVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "start_warm")

make_variable_name(::Type{ColdStartVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "start_cold")

make_variable_name(::Type{EnergyVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "E")

make_variable_name(::Type{EnergyVariableUp}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "Eup")

make_variable_name(::Type{EnergyVariableDown}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "Edown")

make_variable_name(::Type{EnergySurplusVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "energy_surplus")

make_variable_name(::Type{EnergyShortageVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "energy_shortage")

make_variable_name(::Type{LiftVariable}) = :lift

make_variable_name(::Type{LiftVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "lift")

make_variable_name(::Type{OnVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "On")

make_variable_name(::Type{ReactivePowerVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "Q")

make_variable_name(::Type{ReserveVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "R")

make_variable_name(::Type{ServiceRequirementVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "service_requirement")

make_variable_name(::Type{SpillageVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "Sp")

make_variable_name(::Type{StartVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "start")

make_variable_name(::Type{StopVariable}, ::Type{T}) where {T <: PSY.Device} = encode_symbol(T, "stop")

make_variable_name(::Type{SteadyStateFrequencyDeviation}) = :Δf

make_variable_name(::Type{AreaMismatchVariable}) = :area_mismatch

make_variable_name(::Type{DeltaActivePowerUpVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "ΔP_up")

make_variable_name(::Type{DeltaActivePowerUpVariable}, ::Type{PSY.RegulationDevice{T}}) where {T <: PSY.Device} = encode_symbol(T, "ΔP_up")

make_variable_name(::Type{DeltaActivePowerDownVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "ΔP_dn")

make_variable_name(::Type{DeltaActivePowerDownVariable}, ::Type{PSY.RegulationDevice{T}}) where {T <: PSY.Device} = encode_symbol(T, "ΔP_dn")

make_variable_name(::Type{AdditionalDeltaActivePowerUpVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "ΔPe_up")

make_variable_name(::Type{AdditionalDeltaActivePowerUpVariable}, ::Type{PSY.RegulationDevice{T}}) where {T <: PSY.Device} = encode_symbol(T, "ΔPe_up")

make_variable_name(::Type{AdditionalDeltaActivePowerDownVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "ΔPe_dn")

make_variable_name(::Type{AdditionalDeltaActivePowerDownVariable}, ::Type{PSY.RegulationDevice{T}}) where {T <: PSY.Device} = encode_symbol(T, "ΔPe_dn")

make_variable_name(::Type{SmoothACE}, ::Type{T}) where {T <: PSY.AggregationTopology} = encode_symbol(T, "SACE")

make_variable_name(::Type{FlowActivePowerVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "Fp")

make_variable_name(::Type{FlowReactivePowerVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "Fq")

make_variable_name(::Type{FlowActivePowerFromToVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "FpFT")

make_variable_name(::Type{FlowActivePowerToFromVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "FpTF")

make_variable_name(::Type{FlowReactivePowerFromToVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "FqFT")

make_variable_name(::Type{FlowReactivePowerToFromVariable}, ::Type{T}) where {T <: PSY.Component} = encode_symbol(T, "FqTF")


text:
using Potree
using DataStructures
using Visualization
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

function decodeX(x::Int32)::Float64

	potree = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsBIN(v1.7)" # replace this path with local potree directory
	metadata = CloudMetadata(potree)
	boundingBox = metadata.boundingBox
	scale = metadata.scale

	return (x * scale) + boundingBox.x_min
end

function decodeY(x::Int32)::Float64

	potree = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsBIN(v1.7)" # replace this path with local potree directory
	metadata = CloudMetadata(potree)
	boundingBox = metadata.boundingBox
	scale = metadata.scale

	return (x * scale) + boundingBox.y_min
end

function decodeZ(x::Int32)::Float64

	potree = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsBIN(v1.7)" # replace this path with local potree directory
	metadata = CloudMetadata(potree)
	boundingBox = metadata.boundingBox
	scale = metadata.scale

	return (x * scale) + boundingBox.z_min
end

function binToPoints(fname::String)
	data = read(fname)

	#divido in pacchetti di 4 * 4 byte
	treehrc = reshape(data, (16, div(length(data), 16)))

	for i in 1:size(treehrc,2)
		x = parse(Int32, bitstring(UInt8(treehrc[4,i]))*bitstring(UInt8(treehrc[3,i]))*bitstring(UInt8(treehrc[2,i]))*bitstring(UInt8(treehrc[1,i])); base=2)
		y = parse(Int32, bitstring(UInt8(treehrc[8,i]))*bitstring(UInt8(treehrc[7,i]))*bitstring(UInt8(treehrc[6,i]))*bitstring(UInt8(treehrc[5,i])); base=2)
		z = parse(Int32, bitstring(UInt8(treehrc[12,i]))*bitstring(UInt8(treehrc[11,i]))*bitstring(UInt8(treehrc[10,i]))*bitstring(UInt8(treehrc[9,i])); base=2)
		r = parse(Int32, bitstring(UInt8(treehrc[13,i])))
		g = parse(Int32, bitstring(UInt8(treehrc[14,i])))
		b = parse(Int32, bitstring(UInt8(treehrc[15,i])))
		a = parse(Int32, bitstring(UInt8(treehrc[16,i])))

		return vcat(x',y',z'), vcat(r',g',b')

	end
end

function Bpoints(fname::String)
	data = read(fname)

	#divido in pacchetti di 4 * 4 byte
	treehrc = reshape(data, (16, div(length(data), 16)))

	allPoints = Array{Float64,2}(undef, 3, 0)

	for i in 1:size(treehrc,2)
		x = parse(Int32, bitstring(UInt8(treehrc[4,i]))*bitstring(UInt8(treehrc[3,i]))*bitstring(UInt8(treehrc[2,i]))*bitstring(UInt8(treehrc[1,i])); base=2)
		y = parse(Int32, bitstring(UInt8(treehrc[8,i]))*bitstring(UInt8(treehrc[7,i]))*bitstring(UInt8(treehrc[6,i]))*bitstring(UInt8(treehrc[5,i])); base=2)
		z = parse(Int32, bitstring(UInt8(treehrc[12,i]))*bitstring(UInt8(treehrc[11,i]))*bitstring(UInt8(treehrc[10,i]))*bitstring(UInt8(treehrc[9,i])); base=2)

		allPoints = hcat(allPoints, vcat(decodeX(x),decodeY(y),decodeZ(z)))
	end
	return allPoints
end

function Brgb(fname::String)::Lar.Points
	data = read(fname)

	allRGB = Array{LasIO.N0f16,2}(undef, 3, 0)

	#divido in pacchetti di 4 * 4 byte
	treehrc = reshape(data, (16, div(length(data), 16)))

	for i in 1:size(treehrc,2)
		r = parse(UInt32, bitstring(UInt8(treehrc[13,i])))
		g = parse(UInt32, bitstring(UInt8(treehrc[14,i])))
		b = parse(UInt32, bitstring(UInt8(treehrc[15,i])))
		a = parse(UInt32, bitstring(UInt8(treehrc[16,i])))
		allRGB = hcat(allRGB, vcat(r,g,b))

	end
	return allRGB = Array{LasIO.N0f16,2}(undef, 3, 0)
end

# sources
potree = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsBIN(v1.7)" # replace this path with local potree directory
hrcfile = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsBIN(v1.7)/data/r/r.hrc" # hrc

# uso nodi piccoli
binfile = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsBIN(v1.7)/data/r/r0.bin"
lasfile = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsLAS(v1.7)/data/r/r0.las"

metadata = CloudMetadata(potree)
totalPoints = metadata.points
points = Potree.las2larpoints(lasfile)
println(size(points))
rgb = Potree.las2color(lasfile)

pointss = Bpoints(binfile)
rgbs = Brgb(binfile)

println("===========================================")
println("x y z of last point")
println(points[length(points)-2])
println(points[length(points)-1])
println(points[length(points)])
println("===========================================")
println("rgb of last point")
println(rgb[length(rgb)][1])
println("===========================================")

GL.VIEW(
	[
	Visualization.points_color_from_rgb(pointss,rgbs)
	]
)