text:
#=------------------------------------------------------------------------------
                        Formatting Routines
------------------------------------------------------------------------------=#
function produce_ssten_from_triangles(file;use_metis=false,lcc=false)

    A = MatrixNetworks.readSMAT(file)
    (n,m) = size(A)
    if(n != m)
        println("rectangular")
    end

	if !issymmetric(A)
		A = max.(A,A')  #symmetrize for Triangles routine
    end

    if use_lcc
        A, _ = largest_component(A)
    end

	if use_metis
		apply_Metis_permutation!(A)
	end

    T = collect(MatrixNetworks.triangles(A))

	alterfilename = (file,postfix)-> split(file,".smat")[1]*postfix
	if use_metis
	    output_file = alterfilename(file,"_with_metis.ssten")
	else
	    output_file = alterfilename(file,".ssten")
	end

    open(output_file,"w") do f
        write(f,"$(3)\t$(n)\t$(length(T))\n")

        for (v_i,v_j,v_k) in T
            write(f,"$(v_i)\t$(v_j)\t$(v_k)\t1.0\n")
        end
    end

end

function load_UnweightedThirdOrderSymTensor(filepath;enforceFormatting = true)

	#check path validity
	@assert filepath[end-5:end] == ".ssten"

    open(filepath) do file
		#preallocate from the header line
		order, n, m =
			[parse(Int,elem) for elem in split(chomp(readline(file)),'\t')]
        @assert order == 3

		Ti = [ Vector{Tuple{Int,Int}}(undef, 0) for i in 1:n ]

		i = 1
		@inbounds for line in eachline(file)
			entries = split(chomp(line),'\t')

			if enforceFormatting
				(ti,tj,tk) = sort([parse(Int,elem) for elem in entries[1:end-1]])
			else
				(ti,tj,tk) = [parse(Int,elem) for elem in entries[1:end-1]]
			end

			if 0 == ti || 0 == tj || 0 == tk
				error("elements must be indexed by 1.")
			end
			push!(Ti[ti], (tj,tk))
			push!(Ti[tj], (ti,tk))
			push!(Ti[tk], (ti,tj))

		end

		sort!.(Ti)
		return UnweightedThirdOrderSymTensor(n,Ti)
	end

end
"""------------------------------------------------------------------------------
  Loads in a ThirdOrderSymTensor from an ssten file. Data specifications can be
  found in the 'formatting_specification.info' file the data/ folder. The
  enforceFormatting can be used to ensure that the indices are sorted in
  increasing order, and any files which have 0's in indices are updated to be
  indexed by 1.
------------------------------------------------------------------------------"""
function load_ThirdOrderSymTensor(filepath;enforceFormatting = true)

	#check path validity
	@assert filepath[end-5:end] == ".ssten"

    open(filepath) do file
		#preallocate from the header line
		order, n, m =
			[parse(Int,elem) for elem in split(chomp(readline(file)),'\t')]
        @assert order == 3

		indices = Array{Int,2}(undef,order,m)
		values = Array{Float64,1}(undef,m)

		i = 1
		@inbounds for line in eachline(file)
			entries = split(chomp(line),'\t')

			indices[:,i] = [parse(Int,elem) for elem in entries[1:end-1]]
			if enforceFormatting
				sort!(indices[:,i])
			end
			values[i] = parse(Float64,entries[end])
			i += 1
		end

		#check for 0 indexing
		zero_indexed = false

		@inbounds for i in 1:m
		    if indices[1,i] == 0
    			zero_indexed = true
				break
			end
	    end

		if zero_indexed
			indices .+= 1
		end

        return ThirdOrderSymTensor(n,indices,values)
    end
end

function write_to_armadillo_raw_ascii_format(X::Array{T,2},output_file::String) where T
    open(output_file,"w") do f
        for row in eachrow(X)
            println(f,join(row," "))
        end
    end
end

function write_smat(A::SparseMatrixCSC{T,Int},path::String;delimeter::Char=',',kwargs...) where T
	@assert path[end-4:end] == ".smat"
	open(path,"w") do f
		header = join([size(A)...,nnz(A)],delimeter)# ::NTuple{3,Int}
		println(f,header)
		for (i,j,v)=zip(findnz(A)...)
			i -= 1
			j -= 1
			print(f,i,delimeter,j,delimeter)
			println(f,v)
		end
	end
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
module Random

using Base.LibRandom

export srand,
       rand, rand!,
       randn, randn!,
       randbool, randbool!,
       AbstractRNG, RNG, MersenneTwister

abstract AbstractRNG

type MersenneTwister <: AbstractRNG
    state::DSFMT_state
    seed::Union(Uint32,Vector{Uint32})

    function MersenneTwister(seed::Vector{Uint32})
        state = DSFMT_state()
        dsfmt_init_by_array(state, seed)
        return new(state, seed)
    end

    MersenneTwister(seed=0) = MersenneTwister(make_seed(seed))
end

function srand(r::MersenneTwister, seed)
    r.seed = seed
    dsfmt_init_gen_rand(r.state, seed)
    return r
end

## initialization

function __init__()

@unix_only begin
    try
        srand("/dev/urandom")
    catch
        println(STDERR, "Entropy pool not available to seed RNG; using ad-hoc entropy sources.")
        seed = reinterpret(Uint64, time())
        seed = bitmix(seed, uint64(getpid()))
        try
            seed = bitmix(seed, parseint(Uint64, readall(`ifconfig` |> `sha1sum`)[1:40], 16))
        catch
            # ignore
        end
        srand(seed)
    end
end

@windows_only begin
    a = zeros(Uint32, 2)
    win32_SystemFunction036!(a)
    srand(a)
end
end

## srand()

function srand(seed::Vector{Uint32})
    global RANDOM_SEED = seed
    dsfmt_gv_init_by_array(seed)
end
srand(n::Integer) = srand(make_seed(n))

function make_seed(n::Integer)
    n < 0 && throw(DomainError())
    seed = Uint32[]
    while true
        push!(seed, n & 0xffffffff)
        n >>= 32
        if n == 0
            return seed
        end
    end
end

function srand(filename::String, n::Integer)
    open(filename) do io
        a = Array(Uint32, int(n))
        read!(io, a)
        srand(a)
    end
end
srand(filename::String) = srand(filename, 4)

## random floating point values

rand(::Type{Float64}) = dsfmt_gv_genrand_close_open()
rand() = dsfmt_gv_genrand_close_open()

rand(::Type{Float32}) = float32(rand())
rand(::Type{Float16}) = float16(rand())

rand{T<:Real}(::Type{Complex{T}}) = complex(rand(T),rand(T))

rand(r::MersenneTwister) = dsfmt_genrand_close_open(r.state)

## random integers

dsfmt_randui32() = dsfmt_gv_genrand_uint32()
dsfmt_randui64() = uint64(dsfmt_randui32()) | (uint64(dsfmt_randui32())<<32)

rand(::Type{Uint8})   = uint8(rand(Uint32))
rand(::Type{Uint16})  = uint16(rand(Uint32))
rand(::Type{Uint32})  = dsfmt_randui32()
rand(::Type{Uint64})  = dsfmt_randui64()
rand(::Type{Uint128}) = uint128(rand(Uint64))<<64 | rand(Uint64)

rand(::Type{Int8})    = int8(rand(Uint8))
rand(::Type{Int16})   = int16(rand(Uint16))
rand(::Type{Int32})   = int32(rand(Uint32))
rand(::Type{Int64})   = int64(rand(Uint64))
rand(::Type{Int128})  = int128(rand(Uint128))

# Arrays of random numbers

rand(::Type{Float64}, dims::Dims) = rand!(Array(Float64, dims))
rand(::Type{Float64}, dims::Int...) = rand(Float64, dims)

rand(dims::Dims) = rand(Float64, dims)
rand(dims::Int...) = rand(Float64, dims)

rand(r::AbstractRNG, dims::Dims) = rand!(r, Array(Float64, dims))
rand(r::AbstractRNG, dims::Int...) = rand(r, dims)

function rand!{T}(A::Array{T})
    for i=1:length(A)
        A[i] = rand(T)
    end
    A
end
rand(T::Type, dims::Dims) = rand!(Array(T, dims))
rand{T<:Number}(::Type{T}) = error("no random number generator for type $T; try a more specific type")
rand{T<:Number}(::Type{T}, dims::Int...) = rand(T, dims)

## Generate random integer within a range

# remainder function according to Knuth, where rem_knuth(a, 0) = a
rem_knuth(a::Uint, b::Uint) = a % (b + (b == 0)) + a * (b == 0)
rem_knuth{T<:Unsigned}(a::T, b::T) = b != 0 ? a % b : a

# maximum multiple of k <= 2^bits(T) decremented by one,
# that is 0xFFFFFFFF if k = typemax(T) - typemin(T) with intentional underflow
maxmultiple(k::Uint32) = convert(Uint32, div(0x0000000100000000,k + (k == 0))*k - 1)
maxmultiple(k::Uint64) = convert(Uint64, div(0x00000000000000010000000000000000, k + (k == 0))*k - 1)
# maximum multiple of k within 1:typemax(Uint128)
maxmultiple(k::Uint128) = div(typemax(Uint128), k + (k == 0))*k - 1
# maximum multiple of k within 1:2^32 or 1:2^64, depending on size
maxmultiplemix(k::Uint64) = convert(Uint64, div((k >> 32 != 0)*0x0000000000000000FFFFFFFF00000000 + 0x0000000100000000, k + (k == 0))*k - 1)

immutable RandIntGen{T<:Integer, U<:Unsigned}
    a::T   # first element of the range
    k::U   # range length or zero for full range
    u::U   # rejection threshold
end
# generators with 32, 128 bits entropy
RandIntGen{T, U<:Union(Uint32, Uint128)}(a::T, k::U) = RandIntGen{T, U}(a, k, maxmultiple(k))
# mixed 32/64 bits entropy generator
RandIntGen{T}(a::T, k::Uint64) = RandIntGen{T,Uint64}(a, k, maxmultiplemix(k))

# generator for ranges
RandIntGen{T<:Unsigned}(r::UnitRange{T}) = isempty(r) ? error("range must be non-empty") : RandIntGen(first(r), convert(T,last(r) - first(r) + 1))
# specialized versions
for (T, U) in [(Uint8, Uint32), (Uint16, Uint32),
               (Int8, Uint32), (Int16, Uint32), (Int32, Uint32), (Int64, Uint64), (Int128, Uint128),
               (Bool, Uint32), (Char, Uint32)]

    @eval RandIntGen(r::UnitRange{$T}) = isempty(r) ? error("range must be non-empty") : RandIntGen(first(r), convert($U, last(r) - first(r) + 1)) # overflow ok
end

# this function uses 32 bit entropy for small ranges of length <= typemax(Uint32) + 1
# RandIntGen is responsible for providing the right value of k
function rand{T<:Union(Uint64, Int64)}(g::RandIntGen{T,Uint64})
    local x::Uint64
    if (g.k - 1) >> 32 == 0
        x = rand(Uint32)
        while x > g.u
            x = rand(Uint32)
        end
    else
        x = rand(Uint64)
        while x > g.u
            x = rand(Uint64)
        end
    end
    return convert(T, g.a + rem_knuth(x, g.k))
end

function rand{T<:Integer, U<:Unsigned}(g::RandIntGen{T,U})
    x = rand(U)
    while x > g.u
        x = rand(U)
    end
    convert(T, g.a + rem_knuth(x, g.k))
end

rand{T<:Union(Signed,Unsigned,Bool,Char)}(r::UnitRange{T}) = rand(RandIntGen(r))
rand{T<:Real}(r::Range{T}) = convert(T, first(r) + rand(0:(length(r)-1)) * step(r))

function rand!(g::RandIntGen, A::AbstractArray)
    for i = 1 : length(A)
        @inbounds A[i] = rand(g)
    end
    return A
end

rand!{T<:Union(Signed,Unsigned,Bool,Char)}(r::UnitRange{T}, A::AbstractArray) = rand!(RandIntGen(r), A)

function rand!{T<:Real}(r::Range{T}, A::AbstractArray)
    g = RandIntGen(0:(length(r)-1))
    f = first(r)
    s = step(r)
    if s == 1
        for i = 1 : length(A)
            @inbounds A[i] = convert(T, f + rand(g))
        end
    else
        for i = 1 : length(A)
            @inbounds A[i] = convert(T, f + rand(g) * s)
        end
    end
    return A
end

rand{T<:Real}(r::Range{T}, dims::Dims) = rand!(r, Array(T, dims))
rand(r::Range, dims::Int...) = rand(r, dims)

## random Bools

rand!(B::BitArray) = Base.bitarray_rand_fill!(B)

randbool(dims::Dims) = rand!(BitArray(dims))
randbool(dims::Int...) = rand!(BitArray(dims))

randbool() = ((dsfmt_randui32() & 1) == 1)
randbool!(B::BitArray) = rand!(B)

## randn() - Normally distributed random numbers using Ziggurat algorithm

# The Ziggurat Method for generating random variables - Marsaglia and Tsang
# Paper and reference code: http://www.jstatsoft.org/v05/i08/

randn() = randmtzig_randn()
randn(rng::MersenneTwister) = randmtzig_randn(rng.state)
randn!(A::Array{Float64}) = (for i = 1:length(A);A[i] = randmtzig_randn();end;A)
randn!(rng::MersenneTwister, A::Array{Float64}) = (for i = 1:length(A);A[i] = randmtzig_randn(rng.state);end;A)
randn(dims::Dims) = randn!(Array(Float64, dims))
randn(dims::Int...) = randn!(Array(Float64, dims...))

## random UUID generation

immutable UUID
    value::Uint128
end

function uuid4()
    u = rand(Uint128)
    u &= 0xffffffffffff0fff3fffffffffffffff
    u |= 0x00000000000040008000000000000000
    UUID(u)
end

function Base.convert(::Type{Vector{Uint8}}, u::UUID)
    u = u.value
    a = Array(Uint8,36)
    for i = [36:-1:25; 23:-1:20; 18:-1:15; 13:-1:10; 8:-1:1]
        d = u & 0xf
        a[i] = '0'+d+39*(d>9)
        u >>= 4
    end
    a[[24,19,14,9]] = '-'
    return a
end

Base.show(io::IO, u::UUID) = write(io,convert(Vector{Uint8},u))
Base.repr(u::UUID) = ASCIIString(convert(Vector{Uint8},u))

end # module

ext: jl

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
module SkeletonTriangulationsTests

using Test
using Gridap.TensorValues
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Geometry: IN, OUT
using Gridap.Geometry: DiscreteModelMock

model = DiscreteModelMock()

strian = SkeletonTriangulation(model)
test_triangulation(strian)
@test get_background_model(strian) === model
glue = get_glue(strian,Val(1))
@test glue.tface_to_mface === strian.plus.glue.face_to_bgface
glue = get_glue(strian,Val(2))
@test glue.plus.tface_to_mface === strian.plus.glue.face_to_cell
@test glue.minus.tface_to_mface === strian.minus.glue.face_to_cell

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

strian = SkeletonTriangulation(model)
test_triangulation(strian)
@test get_background_model(strian) === model
glue = get_glue(strian,Val(2))
@test glue.tface_to_mface === strian.plus.glue.face_to_bgface
glue = get_glue(strian,Val(3))
@test glue.plus.tface_to_mface === strian.plus.glue.face_to_cell
@test glue.minus.tface_to_mface === strian.minus.glue.face_to_cell

ids = [1,3,4]
vtrian = view(strian,ids)
test_triangulation(vtrian)
@test num_cells(vtrian) == 3
sglue = get_glue(strian,Val(2))
vglue = get_glue(vtrian,Val(2))
@test vglue.tface_to_mface == sglue.tface_to_mface[ids]
sglue = get_glue(strian,Val(3))
vglue = get_glue(vtrian,Val(3))
@test vglue.plus.tface_to_mface == sglue.plus.tface_to_mface[ids]
@test vglue.minus.tface_to_mface == sglue.minus.tface_to_mface[ids]
vn = get_facet_normal(vtrian)
@test isa(vn,SkeletonPair)
@test isa(vn.plus,AbstractArray)
@test isa(vn.minus,AbstractArray)

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model)
Λ = SkeletonTriangulation(Γ)
@test Λ.rtrian === Γ
@test isa(Λ.dtrian,SkeletonTriangulation)
glue = get_glue(Λ,Val(3))
glue = get_glue(Λ,Val(2))
@test is_change_possible(Ω,Λ)
@test is_change_possible(Γ,Λ)

model = DiscreteModelMock()

trian = Triangulation(model)

cell_to_is_left = [true,true,false,true,false]

itrian = InterfaceTriangulation(model,cell_to_is_left)

ltrian = itrian.plus
rtrian = itrian.minus
@test itrian.⁺ === ltrian
@test itrian.⁻ === rtrian

ni = get_facet_normal(itrian)
nl = get_facet_normal(ltrian)
nr = get_facet_normal(rtrian)
@test isa(ni,SkeletonPair)

#using Gridap.Visualization
#
#writevtk(trian,"trian")
#writevtk(itrian,"itrian",nsubcells=10,cellfields=["ni"=>ni,"nl"=>nl,"nr"=>nr])

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

cell_to_inout = zeros(Int8,num_cells(model))
cell_to_inout[1:13] .= OUT
cell_to_inout[14:34] .= IN

itrian = InterfaceTriangulation(model,cell_to_inout)
@test num_cells(itrian) == 11

itrian = InterfaceTriangulation(model,1:13,14:34)
@test num_cells(itrian) == 11
Ω_in = Triangulation(model,findall(i->i==IN,cell_to_inout))
Ω_out = Triangulation(model,findall(i->i==OUT,cell_to_inout))
itrian = InterfaceTriangulation(Ω_in,Ω_out)

#ltrian = get_left_boundary(itrian)
#rtrian = get_right_boundary(itrian)
#ni = get_normal_vector(itrian)
#nl = get_normal_vector(ltrian)
#nr = get_normal_vector(rtrian)
#trian = Triangulation(model)
#
#using Gridap.Visualization
#writevtk(trian,"trian",celldata=["inout"=>cell_to_inout])
#writevtk(itrian,"itrian",nsubcells=10,cellfields=["ni"=>ni,"nl"=>nl,"nr"=>nr])

reffe = LagrangianRefFE(Float64,QUAD,(2,2))
conf = CDConformity((CONT,DISC))
face_own_dofs = get_face_own_dofs(reffe,conf)
strian = SkeletonTriangulation(model,reffe,face_own_dofs)
test_triangulation(strian)
ns = get_facet_normal(strian)
@test length(ns.⁺) == num_cells(strian)

#using Gridap.Visualization
#writevtk(strian,"strian",cellfields=["normal"=>ns])

domain = (0,1,0,1,0,1)
partition = (3,3,3)
oldmodel = CartesianDiscreteModel(domain,partition)
oldstrian = SkeletonTriangulation(oldmodel)

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
Ω1 = Interior(model,[2])
Ω2 = Interior(model,[4])
Γ12 = Interface(Ω1,Ω2)
@test num_cells(Γ12) == 1

end # module


text:
export midpoint, intermediate_point, destination_point, intersection_point,
intersection_points, isinside, isonborder, ison, isselfintersecting,
highest_latitude, highest_latitude_point, lowest_latitude, lowest_latitude_point,
point_cartesian, point_spherical

"""
    midpoint(point₁::Point, point₂::Point)

    output: point::Point

Return the half-way point `midpoint` [deg] on the great circle line between
the positions `point₁` and `point₂` [deg] on a unit sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function midpoint(point₁::Point, point₂::Point)
    Δpos = point₂ - point₁
    Bx = cosd(point₂.ϕ) * cosd(Δpos.λ)
    By = cosd(point₂.ϕ) * sind(Δpos.λ)
    ϕₘ = atand(sind(point₁.ϕ) + sind(point₂.ϕ), √((cosd(point₁.ϕ) + Bx)^2 + By^2))
    λₘ = point₁.λ + atand(By, cosd(point₁.ϕ) + Bx)
    return Point(ϕₘ, normalize(λₘ, -180.0, 180.0))
end

"""
    midpoint(arc::Arc)

    output: point::Point

Return the half-way point `midpoint` [deg] on the great circle `arc`
[deg] on a unit sphere.
"""
midpoint(arc::Arc) = midpoint(arc.point₁, arc.point₂)

"""
    midpoint(arcs::Arcs)

    output: point::Point

Return the half-way point `midpoint` [deg] on the great circle `arcs`
[deg] on a unit sphere.
"""
midpoint(arcs::Arcs) = intermediate_point(arcs, 0.5)

"""
    intermediate_point(point₁::Point, point₂::Point, fraction::Float64)

    output: point::Point

Return the `intermediate_point` [deg] at any `fraction` along the great circle
line section between the two points `point₁` and `point₂` [deg]. The fraction
along the great circle line section is such that `fraction` = 0.0 is at `point₁`
[deg] and `fraction` 1.0 is at `point₂` [deg].

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function intermediate_point(point₁::Point, point₂::Point, fraction::Float64)
    δ = angular_distance(point₁, point₂)
    a = sind((1.0 - fraction) * δ) / sind(δ)
    b = sind(fraction * δ) / sind(δ)
    x = a * cosd(point₁.ϕ) * cosd(point₁.λ) + b * cosd(point₂.ϕ) * cosd(point₂.λ)
    y = a * cosd(point₁.ϕ) * sind(point₁.λ) + b * cosd(point₂.ϕ) * sind(point₂.λ)
    z = a * sind(point₁.ϕ) + b * sind(point₂.ϕ)
    ϕᵢ = atand(z, √(x*x + y*y))
    λᵢ = atand(y, x)
    return Point(ϕᵢ, λᵢ)
end

"""
    intermediate_point(arc::Arc, fraction::Float64)

    output: point::Point

Return the `intermediate_point` [deg] at any `fraction` along the great circle
`arc` [deg]. The fraction along the great circle `arc` is such
that `fraction` = 0.0 is at the start of the `arc` and `fraction` 1.0 is
at the end of the `arc`.
"""
intermediate_point(arc::Arc, fraction::Float64) =
intermediate_point(arc.point₁, arc.point₂, fraction)

"""
    intermediate_point(arcs::Arcs, fraction::Float64)

    output: point::Point

Return the `intermediate_point` [deg] at any `fraction` along the great circle
`arcs` [deg]. The fraction along the great circle `arcs` is such
that `fraction` = 0.0 is at the start of the `arcs` and `fraction` 1.0 is
at the end of the `arcs`.
"""
function intermediate_point(arcs::Arcs, fraction::Float64)
    angular_length_ip = angular_length(arcs) * fraction
    p₁ = arcs.points[1]
    dist = 0.0
    if angular_length_ip == dist
        return arcs.points[1]
    end
    for p₂ in arcs.points[2:end]
        dist_p₁p₂ = angular_distance(p₁, p₂)
        if dist + dist_p₁p₂ ≥ angular_length_ip
            fraction_section = (angular_length_ip - dist) / dist_p₁p₂
            return intermediate_point(p₁, p₂, fraction_section)
        end
        dist += dist_p₁p₂
        p₁ = p₂
    end
    return arcs.points[end]
end

"""
    destination_point(start_point::Point, angular_distance::Float64,
    azimuth::Float64)

    output: point::Point

Given a `start_point` [deg], initial `azimuth` (clockwise from North) [deg],
`angular_distance` [deg] calculate the position of the `destination_point`
[deg] traversing along a great circle line.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function destination_point(start_point::Point, angular_distance::Float64,
azimuth::Float64)
    ϕ₂ = asind(sind(start_point.ϕ) * cosd(angular_distance) +
    cosd(start_point.ϕ) * sind(angular_distance) * cosd(azimuth))
    λ₂ = start_point.λ + atand(sind(azimuth) * sind(angular_distance) *
    cosd(start_point.ϕ), cosd(angular_distance) - sind(start_point.ϕ) * sind(ϕ₂))
    return Point(ϕ₂, λ₂)
end

"""
    destination_point(line::Line, angular_distance::Float64)

    output: point::Point

Given a line with a start_point and azimuth, calculate the `destination_point`
[deg] at the `angular_distance`.
"""
destination_point(line::Line, angular_distance::Float64) =
destination_point(line.point, angular_distance, line.azimuth)

"""
    intersection_point(point₁::Point, point₂::Point, azimuth₁₃::Float64,
    azimuth₂₃::Float64)

    output: point::Point

Return the intersection point `point₃` [deg] of two great circle lines given
by two start points [deg] `point₁` and `point₂` [deg] and two azimuths [deg]
from `point₁` to `point₃` [deg], and from `point₂` to `point₃` [deg].

Under certain circumstances the results can be an ∞ or *ambiguous solution*.

Source: edwilliams.org/avform.htm
"""
function intersection_point(point₁::Point, point₂::Point, azimuth₁₃::Float64,
azimuth₂₃::Float64)
    Δpos = point₂ - point₁
    δ₁₂ = 2asind(√(sind(Δpos.ϕ / 2)^2 + cosd(point₁.ϕ) * cosd(point₂.ϕ) * sind(Δpos.λ/2)^2))
    if sind(δ₁₂) * cosd(point₁.ϕ) == 0.0 || sind(δ₁₂) * cosd(point₂.ϕ) == 0.0
        return Point(NaN, NaN)
    end
    θₐ = acosd((sind(point₂.ϕ) - sind(point₁.ϕ) * cosd(δ₁₂)) / (sind(δ₁₂) * cosd(point₁.ϕ)))
    θᵦ = acosd((sind(point₁.ϕ) - sind(point₂.ϕ) * cosd(δ₁₂)) / (sind(δ₁₂) * cosd(point₂.ϕ)))
    local θ₁₂, θ₂₁
    if sind(Δpos.λ) > 0.0
        θ₁₂ = θₐ
        θ₂₁ = 360.0 - θᵦ
    else
        θ₁₂ = 360.0 - θₐ
        θ₂₁ = θᵦ
    end
    α₁ = (azimuth₁₃ - θ₁₂ + 180.0)%(360.0) - 180.0
    α₂ = (θ₂₁ - azimuth₂₃ + 180.0)%(360.0) - 180.0

    if sind(α₁) ≈ 0.0 && sind(α₂) ≈ 0.0
        return Point(Inf, Inf)  # infinity of intersections
    elseif sind(α₁) * sind(α₂) < 0.0
        # return Point_rad(Inf, Inf)  # intersection ambiguous
        return Point(NaN, NaN)
    else
        α₃ = acosd(-cosd(α₁) * cosd(α₂) + sind(α₁) * sind(α₂) * cosd(δ₁₂))
        δ₁₃ = atand(sind(δ₁₂) * sind(α₁) * sind(α₂), cosd(α₂) + cosd(α₁) * cosd(α₃))
        ϕ₃ = asind(sind(point₁.ϕ) * cosd(δ₁₃) + cosd(point₁.ϕ) * sind(δ₁₃) * cosd(azimuth₁₃))
        Δλ₁₃ = atand(sind(azimuth₁₃) * sind(δ₁₃) * cosd(point₁.ϕ), cosd(δ₁₃)
        - sind(point₁.ϕ) * sind(ϕ₃))
        λ₃ = point₁.λ + Δλ₁₃
        return Point(ϕ₃, λ₃)
    end
end

"""
    intersection_point(line₁::Line, line₂::Line)

    output: point::Point

Return the intersection point `point₃` [deg] of two great circle lines 'line₁'
and `line₂`.

Under certain circumstances the results can be an ∞ or *ambiguous solution*.
"""
intersection_point(line₁::Line, line₂::Line) = intersection_point(line₁.point,
line₂.point, line₁.azimuth, line₂.azimuth)

"""
    intersection_point(point₁::Point, point₂::Point, point₃::Point, point₄::Point)

    output: point::Point

Return the intersection `point` [deg] of two great circle line section given
by two sets of two points: `point₁` and `point₂` [deg] and `point₃` and
`point₄` [deg].

Under certain circumstances the results can be an ∞ or *ambiguous solution*.
https://blog.mbedded.ninja/mathematics/geometry/spherical-geometry/finding-the-intersection-of-two-arcs-that-lie-on-a-sphere/
"""
function intersection_point(point₁::Point, point₂::Point, point₃::Point, point₄::Point)
    if point₁ == point₃ || point₁ == point₄
        return point₁
    elseif point₂ == point₃ || point₂ == point₄
        return point₂
    end

    sph_point₁₁ = point_cartesian(point₁)
    sph_point₁₂ = point_cartesian(point₂)
    sph_point₂₁ = point_cartesian(point₃)
    sph_point₂₂ = point_cartesian(point₄)
    N₁ = cross(sph_point₁₁,sph_point₁₂)
    N₂ = cross(sph_point₂₁, sph_point₂₂)
    L = cross(N₁, N₂)
    I₁ = L/norm(L)
    I₂ = -I₁

    angular_dist(sph_point₁::Vector{Float64}, sph_point₂::Vector{Float64}) =
    angular_distance(point_spherical(sph_point₁), point_spherical(sph_point₂))

    for i in 1:2
        i == 1 ? Iₙ = I₁ : Iₙ = I₂
        θ₁₁ᵢₙ = angular_dist(sph_point₁₁, Iₙ)
        θ₁₂ᵢₙ = angular_dist(sph_point₁₂, Iₙ)
        θ₁₁₁₂ = angular_dist(sph_point₁₁, sph_point₁₂)
        θ₂₁ᵢₙ = angular_dist(sph_point₂₁, Iₙ)
        θ₂₂ᵢₙ = angular_dist(sph_point₂₂, Iₙ)
        θ₂₁₂₂ = angular_dist(sph_point₂₁, sph_point₂₂)
        if abs(θ₁₁ᵢₙ + θ₁₂ᵢₙ - θ₁₁₁₂) < tolerance_deg && abs(θ₂₁ᵢₙ + θ₂₂ᵢₙ - θ₂₁₂₂) < tolerance_deg
            return point_spherical(Iₙ)
        end
    end
    return Point(NaN, NaN)
end

"""
    intersection_point(arc₁::Arc, arc₂::Arc)

    output: point::Point

Return the intersection `point` [deg] of two great circle line sections.

Under certain circumstances the results can be an ∞ or *ambiguous solution*.
"""
intersection_point(arc₁::Arc, arc₂::Arc) =
intersection_point(arc₁.point₁, arc₁.point₂, arc₂.point₁, arc₂.point₂)

"""
    intersection_points(arcs₁::Arcs, arcs₂::Arcs)

    output: points::Vector of Point

Return the intersection points [deg] of two line sections.
"""
function intersection_points(arcs₁::Arcs, arcs₂::Arcs)
    lss₁_p₁ = arcs₁.points[1]
    inter_points = Vector{Point{Float64}}()
    for lss₁_p₂ in arcs₁.points[2:end]
        lss₂_p₁ = arcs₂.points[1]
        distances_section = Vector{Float64}()
        inter_points_section = Vector{Point{Float64}}()
        for lss₂_p₂ in arcs₂.points[2:end]
            intersection_pnt = intersection_point(lss₁_p₁, lss₁_p₂, lss₂_p₁, lss₂_p₂)
            if isinf(intersection_pnt.ϕ) || isnan(intersection_pnt.ϕ)
                #
            elseif intersection_pnt in inter_points_section || intersection_pnt in inter_points
                #
            else
                ipin = false
                for ip in inter_points_section
                    if isapprox(intersection_pnt.ϕ, ip.ϕ; atol=tolerance_deg) &&
                        isapprox(intersection_pnt.λ, ip.λ; atol=tolerance_deg)
                        ipin = true
                        break
                    end
                end
                for ip in inter_points
                    if isapprox(intersection_pnt.ϕ, ip.ϕ; atol=tolerance_deg) &&
                        isapprox(intersection_pnt.λ, ip.λ; atol=tolerance_deg)
                        ipin = true
                        break
                    end
                end
                if ipin
                    #
                else
                    append!(distances_section,
                    angular_distance(lss₁_p₁, intersection_pnt))
                    append!(inter_points_section, [intersection_pnt])
                end
            end
            lss₂_p₁ = lss₂_p₂
        end
        lss₁_p₁ = lss₁_p₂
        p = sortperm(distances_section)
        append!(inter_points, inter_points_section[p])
    end
    return unique(inter_points)
end

intersection_points(arc::Arc, arcs::Arcs) =
intersection_points(Arcs([arc.point₁, arc.point₂]), arcs)

"""
    intersection_points(arcs::Arcs, polygon::Polygon)

    output: points::Vector of Point

Return the intersection points [deg] of a `arcs` with a `polygon`.
"""
intersection_points(arcs::Arcs, polygon::Polygon) =
intersection_points(arcs, Arcs(polygon.points))

intersection_points(arc::Arc, polygon::Polygon) =
intersection_points(Arcs([arc.point₁, arc.point₂]), polygon)

"""
    intersection_points(polygon₁::Polygon, polygon₂::Polygon)

    output: points::Vector of Point

Return the intersection points [deg] of two polygons.
"""
intersection_points(polygon₁::Polygon, polygon₂::Polygon) =
intersection_points(Arcs(polygon₁.points), Arcs(polygon₂.points))

"""
    self_intersection_points(points::Vector{Point{Float64}})

    output: points::Vector of Point

Return the self intersection points [deg] of a set of points.
"""
function self_intersection_points(points::Vector{Point{Float64}})
    if length(points) < 4
        return Vector{Point{Float64}}()
    end
    lss₁_p₁ = points[1]
    inter_points = Vector{Point{Float64}}()
    for (index, lss₁_p₂) in enumerate(points[2:end-1])
        lss₂_p₁ = points[index+2]
        distances_section = Vector{Float64}()
        inter_points_section = Vector{Point{Float64}}()
        for lss₂_p₂ in points[index+3:end]
            intersection_pnt = intersection_point(lss₁_p₁, lss₁_p₂, lss₂_p₁, lss₂_p₂)
            if isinf(intersection_pnt.ϕ) || isnan(intersection_pnt.ϕ)
                #
            else
                append!(distances_section,
                angular_distance(lss₁_p₁, intersection_pnt))
                append!(inter_points_section, [intersection_pnt])
            end
            lss₂_p₁ = lss₂_p₂
        end
        lss₁_p₁ = lss₁_p₂
        p = sortperm(distances_section)
        append!(inter_points, inter_points_section[p])
    end
    return inter_points
end

self_intersection_points(polygon::Polygon) = self_intersection_points(polygon.points)

"""
    isinside(point::Point, polygon::Polygon)

Determine if the `point` is inside the `polygon`. The border is outside the `polygon`.

Source: M. Bevis and J.L. Chatelain, "Locating a point on a spherical surface
relative to a spherical polygon" 1989
"""
function isinside(point::Point, polygon::Polygon)
    return iseven(length(intersection_points(Arc(point, polygon.inside_point), polygon)))
end

"""
    isinside(arc::Arc, polygon::Polygon)

Determine if the `arc` is fully inside the `polygon`. The border is outside the `polygon`.
"""
function isinside(arc::Arc, polygon::Polygon)
    return isinside(arc.point₁, polygon) && length(intersection_points(arc, polygon)) == 0
end

"""
    isinside(arcs::Arcs, polygon::Polygon)

Determine if the `arcs` is fully inside the `polygon`. The border is outside the `polygon`.
"""
function isinside(arcs::Arcs, polygon::Polygon)
    return isinside(arcs.points[1], polygon) && length(intersection_points(arcs, polygon)) == 0
end

"""
    isinside(polygon₁::Polygon, polygon₂::Polygon)

Determine if the `polygon₁` is fully inside the `polygon₂`. The border is outside the `polygon₁`.
"""
function isinside(polygon₁::Polygon, polygon₂::Polygon)
    return isinside(polygon₁.points[1], polygon₂) && length(intersection_points(polygon₁, polygon₂)) == 0
end

"""
    ison(point::Point, line::Line [,tolerance::Float64=tolerance_deg])

Determine if the `point` is on the `line`.
"""
function ison(point::Point, line::Line, tolerance::Float64=tolerance_deg)
    return abs(angular_distance(point, line)) < tolerance
end

"""
    ison(point::Point, arc::Arc [,tolerance::Float64=tolerance_deg])

Determine if the `point` is on the `arc`.
"""
function ison(point::Point, arc::Arc, tolerance::Float64=tolerance_deg)
    return abs(angular_distance(point, arc)) < tolerance
end

"""
    ison(point::Point, arc::Arcs [,tolerance::Float64=tolerance_deg])

Determine if the `point` is on the `arcs`.
"""
function ison(point::Point, arcs::Arcs, tolerance::Float64=tolerance_deg)
    return abs(angular_distance(point, arcs)) < tolerance
end

"""
    ison(point::Point, polygon::Polygon [,tolerance::Float64=tolerance_deg])

Determine if the `point` is on the border of the `polygon`.
"""
function ison(point::Point, polygon::Polygon, tolerance::Float64=tolerance_deg)
    return abs(angular_distance(point, polygon)) < tolerance
end

"""
    along_line_point(point₃::Point, line₁::Line)

    output: point::Point

The `along_line_point` from the start `point₁` [deg] to the closest
point on the great circle line `line₁` to `point₃`.
"""
function along_line_point(point₃::Point, line₁::Line)
    along_dist = along_line_angular_distance(point₃, line₁)
    p3a = destination_point(line₁, along_dist)
    p3b = destination_point(line₁, -along_dist)
    dist_p3a = angular_distance(point₃, p3a)
    dist_p3b = angular_distance(point₃, p3b)
    return dist_p3a < dist_p3b ? p3a : p3b
end

"""
    isselfintersecting(polygon::Polygon)
"""
function isselfintersecting(polygon::Polygon)
    polypoints = deepcopy(polygon.points)
    fraction = 0.999999
    updatelast = intermediate_point(polypoints[end-1], polypoints[end], fraction)
    polypoints[end] = updatelast
    isselfintersecting(polypoints)
end

"""
    isselfintersecting(arcs::Arcs)
"""
isselfintersecting(arcs::Arcs) = isselfintersecting(arcs.points)

"""
    isselfintersecting(points::Vector{Point{Float64}})
"""
function isselfintersecting(points::Vector{Point{Float64}})
    inter_pnts = self_intersection_points(points)
    return length(inter_pnts) != 0
end

"""
    highest_latitude(point::Point, azimuth::Float64)

    output: latitude [deg]

The highest latitude a great circle reaches when starting in a point along the
given azimuth [deg].

Source: https://www.edwilliams.org/avform147.htm#Clairaut
"""
function highest_latitude(point₁::Point, azimuth₁₂::Float64)
    return acosd(abs(sind(azimuth₁₂)*cosd(point₁.ϕ)))
end

"""
    lowest_latitude(point::Point, azimuth::Float64)

    output: latitude [deg]

The lowest latitude a great circle reaches when starting in a point along the
given azimuth [deg].
"""
function lowest_latitude(point₁::Point, azimuth₁₂::Float64)
    return -highest_latitude(point₁, azimuth₁₂)
end

"""
    highest_latitude_point(point::Point, azimuth::Float64)

    output: point::Point

The highest latitude point a great circle reaches when starting in point along the
given azimuth [deg].

Source: https://www.edwilliams.org/avform147.htm#XTE
"""
function highest_latitude_point(point₁::Point, azimuth₁₂::Float64)
    pole = Point(90.0, 0.0)
    Δazimuth = azimuth(point₁,pole) - azimuth₁₂
    dist₁ₚ = angular_distance(point₁, pole)
    p = atand(sind(dist₁ₚ)*cosd(Δazimuth), cosd(dist₁ₚ))
    return destination_point(point₁, p, azimuth₁₂)
end

"""
    lowest_latitude_point(point::Point, azimuth::Float64)

    output: point::Point

The lowest latitude point a great circle reaches when starting in point along the
given azimuth [deg].

Source: https://www.edwilliams.org/avform147.htm#XTE
"""
function lowest_latitude_point(point₁::Point, azimuth₁₂::Float64)
    hlp = highest_latitude_point(point₁, azimuth₁₂)
    return Point(-hlp.ϕ, normalize(hlp.λ+180.0, -180.0, 180.0))
end

#TBD Test scenarios for highest_latitude_point
"""
    highest_latitude_point(point₁::Point, point₂::Point)

    output: point::Point

The highest latitude point a great circle reaches when starting in point₁
towards point₂. The point must be between point₁ and point₂.

Source: https://www.edwilliams.org/avform147.htm#XTE
"""
function highest_latitude_point(point₁::Point, point₂::Point)
    azimuth₁₂ = azimuth(point₁, point₂)
    hlp = highest_latitude_point(point₁, azimuth₁₂)
    if point₁.λ == point₂.λ
        if point₁.ϕ > point₂.ϕ
            return point₁
        else
            return point₂
        end
    else
        if point₁.λ ≤ hlp.λ ≤ point₂.λ || point₂.λ ≤ hlp.λ ≤ point₁.λ
            return hlp
        else
            if point₁.ϕ > point₂.ϕ
                return point₁
            else
                return point₂
            end
        end
    end
end

"""
    lowest_latitude_point(point₁::Point, point₂::Point)

    output: point::Point

The lowest latitude point a great circle reaches when starting in point₁
towards point₂. The point must be between point₁ and point₂.

Source: https://www.edwilliams.org/avform147.htm#XTE
"""
function lowest_latitude_point(point₁::Point, point₂::Point)
    azimuth₁₂ = azimuth(point₁, point₂)
    llp = lowest_latitude_point(point₁, azimuth₁₂)
    if point₁.λ == point₂.λ
        if point₁.ϕ > point₂.ϕ
            return point₂
        else
            return point₁
        end
    else
        if point₁.λ ≤ llp.λ ≤ point₂.λ || point₂.λ ≤ llp.λ ≤ point₁.λ
            return llp
        else
            if point₁.ϕ > point₂.ϕ
                return point₂
            else
                return point₁
            end
        end
    end
end

point_cartesian(point::Point) = [cosd(point.λ)*cosd(point.ϕ); sind(point.λ)*cosd(point.ϕ); sind(point.ϕ)]
point_spherical(point::Vector) = Point(asind(point[3]), atand(point[2], point[1]))


text:
module Milkbot

using  CSV, HTTP, StructTypes,  Plots, DataFrames, JSON3, TensorCast, Setfield, Printf
export MilkPoint, Lactation, FittedLactation, LactationSet, MBParams, Fit, Priors, FitOptions, FittingJob, Credentials  #, StructTypes.StructType
export getPriorsTree, fitLactation, fitLactationSet, mbParamsArray, milkbot, lactationPlot, paramsPlot, refitPlot, fitPlot

struct Credentials
    server::String
    apikey::String
end

struct MilkPoint
    dim::Integer
    milk::Float64
end
MilkPoint((dim,milk)) = MilkPoint(dim,milk)

struct Lactation
    lacKey::String
    breed::String
    parity::Integer
    points
end

struct MBParams
    scale::Float16
    ramp::Float16
    decay::Float16
    offset::Float16
end
MBParams(js ::JSON3.Object{Base.CodeUnits{UInt8, String}, SubArray{UInt64, 1, Vector{UInt64}, Tuple{UnitRange{Int64}}, true}}) = MBParams(js.scale,js.ramp,js.decay,js.offset)
MBParams(d::Dict) = MBParams(d["scale"],d["ramp"],d["decay"],d["offset"])

const mbParamsArray(p::MBParams) = [p.scale, p.ramp,p.decay,p.offset]
struct Fit
    lacKey::String
    params::MBParams
    n
    sdResiduals
    path
end
struct FittedLactation
    lac:: Lactation
    fit:: Fit
end
struct Priors
    means::MBParams
    sd::MBParams
    seMilk::Float16
    milkUnit::String
end
struct FitOptions
    steppedFit::Bool
    returnInputData::Bool
    returnPriors::Bool
    returnPath::Bool
    fitMethod::String  #: String = "MilkBot@2.01",
    preferredMilkUnit::String #: MilkUnit = Kg
end
const trackedOptions = FitOptions(false,false,false,true,"MilkBot@2.01","Kg")
struct LactationSet
    name #: Option[String],
    lactations #: Seq[Lactation],
    milkUnit #: Option[MilkUnit.Value]
end
struct FittingJob
    lactationSet::LactationSet
    priorsTree
    options
end
const TrackedFittingJob(l::Lactation) = FittingJob(LactationSet([l]),nothing,trackedOptions)

StructTypes.StructType(::Type{FittingJob}) = StructTypes.Struct()
StructTypes.StructType(::Type{LactationSet}) = StructTypes.Struct()
StructTypes.StructType(::Type{Lactation}) = StructTypes.Struct()
StructTypes.StructType(::Type{MilkPoint}) = StructTypes.Struct()
StructTypes.StructType(::Type{FitOptions}) = StructTypes.Struct()
StructTypes.StructType(::Type{Priors}) = StructTypes.Struct()
StructTypes.StructType(::Type{MBParams}) = StructTypes.Struct()
StructTypes.StructType(::Type{Fit}) = StructTypes.Struct()

const getPriorsTree(cred::Credentials) =
begin
    try
        response = HTTP.get("$(cred.server)/priorstree?milkUnit=Kg", Dict("X-API-KEY"=>cred.apikey))
        return String(response.body)
    catch e
        return "Error occurred : $e"
    end
end

const fitLactation(lac::Lactation, cred::Credentials) =
begin
    jsbody = JSON3.write("lactation" => lac)
    try
        response = HTTP.post("$(cred.server)/fitlactation", Dict("X-API-KEY"=>cred.apikey),jsbody)
        return String(response.body)
    catch e
        return "Error occurred : $e"
    end
end

const fitLactation(lac::Lactation, priors::Priors, cred::Credentials) =
begin
    jsbody = JSON3.write(Dict("lactation" => lac, "priors" => priors))
    try
        response = HTTP.post("$(cred.server)/fitlactation", Dict("X-API-KEY"=>cred.apikey),jsbody)
        return String(response.body)
    catch e
        return "Error occurred : $e"
    end
end

const fitLactation(lac::Lactation, priors::Priors, cred::Credentials, withPath::Bool) =
begin
    jsbody = JSON3.write(Dict("lactation" => lac, "priors" => priors))
    try
        response = HTTP.post("$(cred.server)/fitlactation?includePath=$(withPath)", Dict("X-API-KEY"=>cred.apikey),jsbody)
        return String(response.body)
    catch e
        return "Error occurred : $e"
    end
end

const fitLactationSet(fj::FittingJob, cred::Credentials) =
begin
    jsbody = JSON3.write("fittingJob" => fj)
    try
        response = HTTP.post("$(cred.server)/fitlactations", Dict("X-API-KEY"=>cred.apikey),jsbody)
        return String(response.body)
    catch e
        return "Error occurred : $e"
    end
end

#const milkbot(r::RowTypeSitko) = t ->  ((1 - ℯ^((r.offset - t)/r.ramp)/2)*r.scale)/ℯ^(r.decay * t)
const milkbot(f::Fit) = t ->  ((1 - ℯ^((f.params.offset - t)/f.params.ramp)/2)*f.params.scale)/ℯ^(f.params.decay * t)
const milkbot(p:: MBParams) = t ->  ((1 - ℯ^((p.offset - t)/p.ramp)/2)*p.scale)/ℯ^(p.decay * t)

const lactationPlot(fl::FittedLactation) =
 begin
    dim = map(m -> m.dim, fl.lac.points)
    milk = map(m -> m.milk, fl.lac.points)
    plot(dim  ,milk, seriestype = :scatter, color = :black, label = fl.lac.lacKey)
    plot!(milkbot(fl.fit),0,maximum(dim), label = "fit")
 end

 const paramsPlot(fit::Fit, priors::Priors) =
 begin
     data = [
        map(d->(d["scale"]-priors.means.scale)/priors.sd.scale,fit.path),
        map(d->(d["ramp"]-priors.means.ramp +.004)/priors.sd.ramp,fit.path),
        map(d->(d["decay"]-priors.means.decay)/priors.sd.decay,fit.path),
        map(d->(d["offset"]-priors.means.offset-.004)/priors.sd.offset-.01,fit.path)
     ]
      plot(1:length(fit.path),data , labels = ["scale" "ramp" "decay" "offset"], legend=:topleft)
 end

 fitPlot(l::Lactation, priors::Priors, cred:: Credentials) =
 begin
    dim = map(m -> m.dim, l.points)
    milk = map(m -> m.milk, l.points)
    fit = JSON3.read(fitLactation(l, priors, cred, true), Fit) #call to API to re-fit
    plot(dim,milk, seriestype = :scatter, label = l.lacKey)
    s = map(p ->milkbot(MBParams(p)) , collect(fit.path)) #vector of mb functions
    plot!(pop!(s),0,maximum(dim), color = :black, linewidth =2, label = "final",
    annotate = (((.4, .15), (@sprintf("SE: %.2f", fit.sdResiduals), 10,:right))))
    p1 = plot!(s[1:min(12, length(s))],0,maximum(dim), label = "")
    p2 =   paramsPlot(fit, priors)
    df = vcat(DataFrame.(fit.path)...)
    plot(p1,p2, layout = (2,1), size = (600,600))
    #vcat(DataFrame.(reFit.path)...)
 end

 #version including API call
 const refitPlot(fl::FittedLactation, priors::Priors, cred:: Credentials) =
 begin
    dim = map(m -> m.dim, fl.lac.points)
    milk = map(m -> m.milk, fl.lac.points)
    reFit = JSON3.read(fitLactation(fl.lac, priors, cred, true), Fit) #call to API to re-fit
    plot(dim, milk, seriestype = :scatter, label = fl.lac.lacKey)
    plot!(milkbot(fl.fit),0,maximum(dim), label = "original fit", linestyle = :dash, color = :black, linewidth =2,
    annotate = (((.4, .25), ("original SE: $(fl.fit.sdResiduals)", 10, :right))))
    s = map(p ->milkbot(MBParams(p)) , collect(reFit.path)) #vector of mb functions
    plot!(pop!(s),0,maximum(dim), color = :black, linewidth =2, label = "re-fit",
    annotate = (((.4, .15), (@sprintf("refit SE: %.2f", reFit.sdResiduals), 10,:right))))
    p1 = plot!(s[1:min(12, length(s))],0,maximum(dim), label = "")
    p2 =   paramsPlot(reFit, priors)
    df = vcat(DataFrame.(reFit.path)...)
    plot(p1,p2, layout = (2,1), size = (600,600))
    #vcat(DataFrame.(reFit.path)...)
 end

  #no API call
 const refitPlot(fl::FittedLactation, reFit:: Fit) =
 begin
    dim = map(m -> m.dim, fl.lac.points)
    milk = map(m -> m.milk, fl.lac.points)
    plot(dim  , milk, seriestype = :scatter, label = fl.lac.lacKey)
    plot!(milkbot(fl.fit),0,maximum(dim), label = "original fit", linestyle = :dash, color = :black, linewidth =2)
    s = map(p ->milkbot(MBParams(p)) , collect(reFit.path)) #vector of mb functions
    plot!(pop!(s),0,maximum(dim), color = :black, linewidth =2, label = "re-fit")
    p1 = plot!(s[1:min(12, length(s))],0,maximum(dim), label = "")
    p2 =   paramsPlot(reFit, priors)
    df = vcat(DataFrame.(reFit.path)...)
    plot(p1,p2, layout = (2,1), size = (600,600))
 end

end #Milkbot module

max_stars_repo_name: jehrlich/MilkBot-API

max_issues_repo_name: jehrlich/MilkBot-API

max_forks_repo_name: jehrlich/MilkBot-API


text:
module Testing
  using CxxWrap, MPI
  import ..libjltrilinos
  import ..Teuchos

  @wrapmodule(libjltrilinos, :register_testing)
end


text:
using Statistics
using Distributions
using ProgressMeter

# Defining the abstract node type
abstract type node end

# Defining structures
mutable struct tree_node <: node
        node_label               :: Int #Node Label
        observations_index       :: Array{Int64,1} # Index of Observations
        depth                    :: Int64 # Depth of the Tree
        node_var                 :: Int64 # Index of Splitted var
        node_var_split           :: Float64 # Split var (first we will considerate continuous case)
        terminal                 :: Bool  # If is terminal or not
        parent_node              :: Int64 # Index of parent node
        mu                       :: Float64 # Mu value of distribution of observations that compose that node
        left_node                :: Bool
        right_node               :: Bool
end

mutable struct root <: node
        x                        :: AbstractMatrix # Index of Observations
        depth                    :: Int64 # Depth of the Tree
        mu                       :: Float64 # Mu value of distribution of observations that compose that node
end

mutable struct tree
        root    :: node
        nodes   :: Array{node,1}
end

#Changing print messages

function Base.show(io::IO,r::root)
    print(io,"Root node with ",size(r.x)[1]," observations.\n",
    " μ: ",r.mu,"\n")
end

function Base.show(io::IO,t::tree)
    print(io,"Tree with ",length(t.nodes)," nodes")
end

# Think if need to print more observations
function Base.show(io::IO,n::tree_node)
    print(io,"Node Label: ",n.node_label,"\n",
            " Observations: ",n.observations_index,"\n",
            " Depth: ",n.depth,"\n",
            " Node Var: ",n.node_var,"\n",
            " Node Var Split: ",n.node_var_split,"\n",
            " Parent Node: ",n.parent_node,"\n",
            " Is terminal: ",n.terminal,"\n",
            " μ: ",n.mu,"\n"
    )
end

#==== BART FUNCTIONS ====#

# Grow for root
function grow(node_to_be_spliten::root,min_node_size)

    # Creating a mecanism to reject trees that produce nodes with n<min_node_size
    bad_tree=true
    count_bad_tree=0

    # Creating the left and right node
    left_node_index=Int64[]
    right_node_index=Int64[]
    p_split = Int64
    p_split_value= Float64

    while bad_tree

        #Cleaning the nodes
        left_node_index=[]
        right_node_index=[]

        p_split=rand(1:size(node_to_be_spliten.x)[2],1)[1] # Choose what variable will
                                                                              # be used to split the root node
        p_split_value=rand(node_to_be_spliten.x[:,p_split],1)[1] # Selecting the splitting rule for that var

        # Splitting the node
        for i in 1:size(node_to_be_spliten.x)[1]
            if node_to_be_spliten.x[i,p_split]<p_split_value
                push!(left_node_index,i)
            else
                push!(right_node_index,i)
            end
        end

        #Verying the length of the new terminal nodes
        if (length(left_node_index)<min_node_size) || (length(right_node_index)<min_node_size)
            count_bad_tree+=1
        else
            bad_tree=false
        end

        # To not be trapped in a loop
        if count_bad_tree==2
            return tree(node_to_be_spliten,tree_node[])
        end

    end

    # Defining the children
    left_child_node=tree_node(1,left_node_index,1,p_split,p_split_value,true,0,0,true,false)
    right_child_node=tree_node(2,right_node_index,1,p_split,p_split_value,true,0,0,false,true)

    return tree(node_to_be_spliten,node[left_child_node,right_child_node])

end

# Grow for tree
function grow(tree_to_be_grown::tree,min_node_size)

    # Creating the new tree
    new_tree=deepcopy(tree_to_be_grown)

    terminal_nodes=node[]
    terminal_node_label=Int
    # Sweep all tree nodes
    for i in new_tree.nodes
        if i.terminal #Verify if it is a terminal node
            push!(terminal_nodes,i) # Add the terminal node
        end
    end

    terminal_to_be_spliten_index=rand(1:length(terminal_nodes),1)

    terminal_node_to_be_spliten=terminal_nodes[terminal_to_be_spliten_index][1]

    # Creating a mecanism to reject trees that produce nodes with n<min_node_size
    bad_tree=true
    count_bad_tree=0

    # Creating the left and right node
    left_node_index=Int64[]
    right_node_index=Int64[]
    p_split = Int64
    p_split_value= Float64

    while bad_tree

        #Cleaning the nodes
        left_node_index=[]
        right_node_index=[]

        p_split=rand(1:size(new_tree.root.x)[2],1)[1] # Choose what variable will
                                                                              # be used to split the root node
        # Selecting the splitting rule for that var
        p_split_value=rand(new_tree.root.x[terminal_node_to_be_spliten.observations_index,p_split],1)[1]

        # Splitting the node
        for i in terminal_node_to_be_spliten.observations_index
            if new_tree.root.x[i,p_split]<p_split_value
                push!(left_node_index,i)
            else
                push!(right_node_index,i)
            end
        end

        #Verying the length of the new terminal nodes
        if (length(left_node_index)<min_node_size) || (length(right_node_index)<min_node_size)
            count_bad_tree+=1
        else
            bad_tree=false
        end

        # To not be trapped in a loop
        if count_bad_tree==2
            return tree_to_be_grown
        end

    end

    # Defining the children
    left_child_node=tree_node(length(tree_to_be_grown.nodes)+1,
                                left_node_index,terminal_node_to_be_spliten.depth+1,p_split,
                                p_split_value,
                                true,terminal_node_to_be_spliten.node_label,0,true,false)
    right_child_node=tree_node(length(tree_to_be_grown.nodes)+2,
                                right_node_index,terminal_node_to_be_spliten.depth+1,p_split,
                                p_split_value,
                                true,terminal_node_to_be_spliten.node_label,0,false,true)

    # Updating the splited node
     for i in new_tree.nodes
        if(i.node_label==terminal_node_to_be_spliten.node_label)
            i.terminal=false
            break
        end
     end

    # Adding the new nodes
    push!(new_tree.nodes,left_child_node,right_child_node)

    return new_tree

end

# Prune

function prune(tree_to_be_pruned::tree)

     # Loading the new tree
    new_tree=deepcopy(tree_to_be_pruned)

    # Case of a root node only
    if(length(new_tree.nodes)==0)
        return new_tree
    end

    # In case of a shallow tree (just two terminal nodes)
    if(length(new_tree.nodes)==2)
        return tree(new_tree.root,tree_node[])
    end

    parent_nodes_indexes=Int[]
    parent_to_be_pruned_index=Int
    child_node_pruned=Int
    children_index=Int[0,0]

    # Sweep all tree nodes
    for i in 1:length(new_tree.nodes)
        if !new_tree.nodes[i].terminal #Verify if it is not a terminal node
            push!(parent_nodes_indexes,i) # Add the parent_node
        end
    end

    bad_to_prune=true
    # Veryfing if its children are both terminal nodes

    while(bad_to_prune)
        terminal_counter=0
        children_index=[0,0]

        parent_to_be_pruned_index=rand(parent_nodes_indexes,1)[1]

        children=get_children(new_tree,new_tree.nodes[parent_to_be_pruned_index])

        # Here I verify if the parent node is the chosen one, and if the node is terminal
        if(length(children)==2)
            children_index[1]=children[1]
            children_index[2]=children[2]
            bad_to_prune=false
        end

    end

    # Transforming back as terminal
    new_tree.nodes[parent_to_be_pruned_index].terminal=true

    # Removing the children

    deleteat!(new_tree.nodes,children_index)

    return new_tree

end

# Change function

function change(tree_to_be_changed::tree,min_node_size)
    new_tree=deepcopy(tree_to_be_changed)

    # If the tree is just composed by its terminal node

    if(length(tree_to_be_changed.nodes)==0)
        return new_tree
    end

    parent_nodes_indexes=Int[]
    parent_to_be_pruned_index=Int
    child_node_pruned=Int
    children_index=Int[]
    # Sweep all tree nodes
    for i in 1:length(new_tree.nodes)
        if !new_tree.nodes[i].terminal #Verify if it is not a terminal node
            push!(parent_nodes_indexes,i) # Add the parent_node
        end
    end

    # Creating a mecanism to reject trees that produce nodes with n<min_node_size
    bad_tree=true
    count_bad_tree=0
    new_tree_updated_list=Any[]

    while(bad_tree)

        # Selecting the node to be changed
        node_to_be_changed_label=rand(push!(parent_nodes_indexes,0),1)[1]

        #print("Node changed: ",node_to_be_changed_label,'\n')

        if(node_to_be_changed_label==0)
            current_node= new_tree.root
        else
            current_node=new_tree.nodes[node_to_be_changed_label]
        end

        # Choose what variable will be used to split the root node
        p_split=rand(1:size(new_tree.root.x)[2],1)[1]

        # Selecting the splitting rule for that var
        if(node_to_be_changed_label==0)
            p_split_value=rand(new_tree.root.x[:,p_split],1)[1]
        else
            p_split_value=rand(new_tree.root.x[current_node.observations_index,p_split],1)[1]
        end

        # Get the chidren from node that will be changed
        children_from_changed_node=get_children(new_tree, current_node)

        # Change the node and split rules
        new_tree.nodes[  children_from_changed_node[1]  ].node_var=p_split # Changing the node var
        new_tree.nodes[  children_from_changed_node[2]  ].node_var=p_split # Changing the node var

        new_tree.nodes[  children_from_changed_node[1]  ].node_var_split =p_split_value  # Changing the node split
        new_tree.nodes[  children_from_changed_node[2]  ].node_var_split =p_split_value #Changing the node split

        # Update the tree with the new splitting rule
        new_tree_updated_list=update_tree_function(new_tree,children_from_changed_node,min_node_size)

        # Verify if its a bad tree
        if(new_tree_updated_list[2])
            count_bad_tree+=1
            new_tree=deepcopy(tree_to_be_changed)
        else
            # Setting bad tree as false therefore go out of the loop
            bad_tree=new_tree_updated_list[2]
        end

        # Limit of tries
        if(count_bad_tree==2)
                return new_tree
        end

    end

    return new_tree_updated_list[1]

end

# Swap

function swap(tree_to_be_swaped::tree,min_node_size)

    new_tree=deepcopy(tree_to_be_swaped)

    parent_nodes_indexes=Int[]
    parent_to_be_pruned_index=Int
    child_node_pruned=Int
    children_index=Int[]

    # Count non-terminal
    for i in 1:length(new_tree.nodes)
        if (!new_tree.nodes[i].terminal && new_tree.nodes[i].parent_node!=0) #Verify if is not terminal child node, where the parent is not the root
            if(!new_tree.nodes[(new_tree.nodes[i].parent_node)].terminal)#Verify if its father is also a terminal node
                push!(parent_nodes_indexes,i) # Add the parent_node
            end
        end
    end

    parent_nodes_indexes=unique(parent_nodes_indexes)

    # If there is no terminal nodes enough
    if(length(parent_nodes_indexes)==0)
        return new_tree
    end

    #print(parent_nodes_indexes)

    # Creating a mecanism to reject trees that produce nodes with n<min_node_size
    bad_tree=true
    count_bad_tree=0
    new_tree_updated_list=Any[]

    while(bad_tree)

        # Selecting the node to be changed
        node_to_be_swaped_label=rand(parent_nodes_indexes,1)[1]

        #print("Node changed: ",node_to_be_swaped_label,'\n')

        # Getting the both internal nodes
        current_node_child=new_tree.nodes[node_to_be_swaped_label]
        current_node_father=new_tree.nodes[current_node_child.parent_node]

        #print(current_node_child,'\n')
        #print(current_node_father,'\n')

        # Get the chidren from node that will be swaped
        children_from_swaped_node_child=get_children(new_tree, current_node_child)[[1,2]]
        children_from_swaped_node_father=get_children(new_tree, current_node_father)

        #print(children_from_swaped_node_child,children_from_swaped_node_father)
        # Swapping the nodes

        new_tree.nodes[ children_from_swaped_node_father[1]].node_var, new_tree.nodes[children_from_swaped_node_child[1]].node_var= new_tree.nodes[children_from_swaped_node_child[1]].node_var,new_tree.nodes[ children_from_swaped_node_father[1]].node_var
        new_tree.nodes[ children_from_swaped_node_father[1]].node_var_split, new_tree.nodes[children_from_swaped_node_child[1]].node_var_split= new_tree.nodes[children_from_swaped_node_child[1]].node_var_split,new_tree.nodes[ children_from_swaped_node_father[1]].node_var_split

        new_tree.nodes[ children_from_swaped_node_father[2]].node_var, new_tree.nodes[children_from_swaped_node_child[2]].node_var= new_tree.nodes[children_from_swaped_node_child[2]].node_var,new_tree.nodes[ children_from_swaped_node_father[2]].node_var
        new_tree.nodes[ children_from_swaped_node_father[2]].node_var_split, new_tree.nodes[children_from_swaped_node_child[2]].node_var_split= new_tree.nodes[children_from_swaped_node_child[2]].node_var_split,new_tree.nodes[ children_from_swaped_node_father[2]].node_var_split

        # Update the tree with the new splitting rule
        new_tree_updated_list=update_tree_function(new_tree,children_from_swaped_node_father,min_node_size)

        #print(new_tree_updated_list[1].nodes)
        # Verify if its a bad tree
        if(new_tree_updated_list[2])
            count_bad_tree+=1
            new_tree=deepcopy(tree_to_be_swaped)
        else
            # Setting bad tree as false therefore go out of the loop
            bad_tree=new_tree_updated_list[2]
        end

        # Limit of tries
        if(count_bad_tree==2)
            #print("BAAAD TREEE")
                return tree_to_be_swaped
        end

    end

    return new_tree_updated_list[1]

end

# Auxiliar functions

# Get the children from a tree node
function get_children(update_tree::tree,parent_node::tree_node,all_children_node_copy=[])

    # Count aux
    count_aux=0
    all_children_node=[]
    for i in 1:length(update_tree.nodes)
    #all_chidren_node=node[]
        if(update_tree.nodes[i].parent_node==parent_node.node_label)
            push!(all_children_node,i)
            #print("Node number ",i,"\n")
            count_aux+=1
        end
    end

    all_children_node_copy=deepcopy(all_children_node)

        if(!update_tree.nodes[all_children_node[end]].terminal)
            push!(all_children_node_copy,get_children(update_tree,update_tree.nodes[all_children_node[end]],all_children_node_copy))
        end

        if(!update_tree.nodes[all_children_node[end-1]].terminal)
            push!(all_children_node_copy,get_children(update_tree,update_tree.nodes[all_children_node[end-1]],all_children_node_copy))
        end

    return collect(Iterators.flatten(all_children_node_copy))

end

# Get children from a root node

function get_children(update_tree::tree,parent_node::root,all_children_node_copy=[])

    # Count aux
    count_aux=0
    all_children_node=[]

    for i in 1:length(update_tree.nodes)
        if(update_tree.nodes[i].parent_node==0)
            push!(all_children_node,i)
            count_aux+=1
        end
    end

    all_children_node_copy=deepcopy(all_children_node)

        if(!update_tree.nodes[all_children_node[end]].terminal)
            push!(all_children_node_copy,get_children(update_tree,update_tree.nodes[all_children_node[end]],all_children_node_copy))
        end

        if(!update_tree.nodes[all_children_node[end-1]].terminal)
            push!(all_children_node_copy,get_children(update_tree,update_tree.nodes[all_children_node[end-1]],all_children_node_copy))
        end

    return collect(Iterators.flatten(all_children_node_copy))

end

# Update a tree

function update_tree_function(tree_to_be_updated::tree,list_of_nodes::Array{Int64,1},min_node_size)

    # Creating a copy of tree that will be updated
    new_tree=deepcopy(tree_to_be_updated)

    # Creating a boolean to indicate if will produce a 'bad_tree' (i.e: few number of observations in terminal node)
    is_a_bad_tree=false

    for i in 1:length(list_of_nodes)

        #print(" Updating node: ",list_of_nodes[i],"\n")
        # Creating the aux of each index
        left_index_aux=Int64[]
        right_index_aux=Int64[]

        if(new_tree.nodes[list_of_nodes[i]].left_node)
            # Consider the case of the parent node is the root node
            if(new_tree.nodes[list_of_nodes[i]].parent_node==0)
                for k in 1:size(new_tree.root.x)[1]
                    if(new_tree.root.x[k,new_tree.nodes[list_of_nodes[i]].node_var]<new_tree.nodes[list_of_nodes[i]].node_var_split)
                        push!(left_index_aux,k)
                    end
                end

                # Left index
                new_tree.nodes[list_of_nodes[i]].observations_index=left_index_aux

                # Validating the tree
                if(length(left_index_aux)<min_node_size)
                    is_a_bad_tree=true
                end

            else
            # Consider the other nodes
            for k in new_tree.nodes[new_tree.nodes[list_of_nodes[i]].parent_node].observations_index
                    if(new_tree.root.x[k,new_tree.nodes[list_of_nodes[i]].node_var]<new_tree.nodes[list_of_nodes[i]].node_var_split)
                        push!(left_index_aux,k)
                    end
            end

            # Left-index
            new_tree.nodes[list_of_nodes[i]].observations_index=left_index_aux

            #print("Left index",left_index_aux,"\n")

            # Validating the tree
            if(length(left_index_aux)<min_node_size)
                is_a_bad_tree=true
            end

            end
        else
            # Consider the case of the parent node is the root node
            if(new_tree.nodes[list_of_nodes[i]].parent_node==0)
                for k in 1:size(new_tree.root.x)[1]
                    if(new_tree.root.x[k,new_tree.nodes[list_of_nodes[i]].node_var]>=new_tree.nodes[list_of_nodes[i]].node_var_split)
                        push!(right_index_aux,k)
                    end
                end

                # Right index
                new_tree.nodes[list_of_nodes[i]].observations_index=right_index_aux

                #print("Right index",left_index_aux,"\n")

                # Validating the tree
                if(length(right_index_aux)<min_node_size)
                    is_a_bad_tree=true
                end
            else
            # Consider the other nodes
                for k in new_tree.nodes[new_tree.nodes[list_of_nodes[i]].parent_node].observations_index
                        if(new_tree.root.x[k,new_tree.nodes[list_of_nodes[i]].node_var]>=new_tree.nodes[list_of_nodes[i]].node_var_split)
                            push!(right_index_aux,k)
                        end
                end

                new_tree.nodes[list_of_nodes[i]].observations_index=right_index_aux

                #print("Right index",left_index_aux,"\n")

                if(length(right_index_aux)<min_node_size)
                        is_a_bad_tree=true
                end

            end
        end

    end
    return [new_tree,is_a_bad_tree]
end

# Moving a tree
function move_tree(tree_to_be_moved::tree,move,min_node_size)

    new_tree=deepcopy(tree_to_be_moved)

    if(move=="grow")

        # In case of a tree with just the root node
        if(length(tree_to_be_moved.nodes)==0)
            new_tree=grow(new_tree.root,min_node_size)
        else
            new_tree=grow(new_tree,min_node_size)
        end

    elseif(move=="prune")
        new_tree=prune(tree_to_be_moved)

    elseif(move=="change")
        new_tree=change(tree_to_be_moved,min_node_size)

    else
        new_tree=swap(tree_to_be_moved,min_node_size)

    end

    return [new_tree,move]

end

# Calculating full conditional
function tree_full_conditional(conditional_tree::tree,residuals::Array{Float64,1},tau::Number,tau_mu::Float64)

    # Consider the case of just the single root node
    if(length(conditional_tree.nodes)==0)
       #Seleting the terminal nodes
           sum_residuals_node_squared=sum(residuals.^2)
           sum_residuals_node=sum(residuals)
           terminal_nodes_size=size(conditional_tree.root.x)[1]

           # See the Math RBART .pdf Equation 1
           log_posterior=0.5*length(residuals)*log(tau) + 0.5*(sum(log.(tau_mu ./(tau_mu.+terminal_nodes_size.*tau))))-
                         0.5*tau*sum(sum_residuals_node_squared)+0.5*(tau^2)*sum((sum_residuals_node.^2) ./(tau_mu.+terminal_nodes_size.*tau))

           return log_posterior

    else

        # Seleting the terminal nodes
        terminal_nodes=tree_node[]
        sum_residuals_node_squared=Float64[]
        sum_residuals_node=Float64[]
        terminal_nodes_size=Int[]

        for i in conditional_tree.nodes
            if(i.terminal)
                push!(terminal_nodes,i)
                push!(sum_residuals_node_squared,sum(residuals[i.observations_index].^2))
                push!(sum_residuals_node,sum(residuals[i.observations_index]))
                push!(terminal_nodes_size,length(i.observations_index))
            end
        end

        # Obtatining the number of terminal nodes
        n_terminal_node=length(terminal_nodes)

        # See the Math RBART .pdf Equation 1
        log_posterior=0.5*length(residuals)*log(tau) + 0.5*(sum(log.(tau_mu ./(tau_mu.+terminal_nodes_size.*tau))))-
                      0.5*tau*sum(sum_residuals_node_squared)+0.5*(tau^2)*sum((sum_residuals_node.^2) ./(tau_mu.+terminal_nodes_size.*tau))

        return log_posterior

    end

end

# Updating the mu values

function update_mu(tree_object::tree,residuals,tau,tau_mu)

    # Creating the dummy tree
    new_tree=deepcopy(tree_object)

    # Case of just have the root node
    if(length(new_tree.nodes)==0)
        mean_values=(tau*sum(residuals))/(size(new_tree.root.x)[1]*tau +tau_mu)
        sd_values=sqrt(1/(size(new_tree.root.x)[1]*tau +tau_mu))
        mu_value=randn(1)[1]*sd_values+mean_values
        new_tree.root.mu=mu_value
        return new_tree

    else

         # Seleting the terminal nodes
        terminal_nodes_index=Int[]
        sum_residuals_node=Float64[]
        terminal_nodes_size=Int[]

        for i in 1:length(new_tree.nodes)
            if(new_tree.nodes[i].terminal)
                push!(terminal_nodes_index,i)
                push!(sum_residuals_node,sum(residuals[new_tree.nodes[i].observations_index]))
                push!(terminal_nodes_size,length(new_tree.nodes[i].observations_index))
            end
        end

        # Obtatining the number of terminal nodes
        n_terminal_node=length(terminal_nodes_index)

        # See the Math RBART .pdf Equation 1
        mean_values=(tau.*sum_residuals_node) ./(terminal_nodes_size*tau .+tau_mu)
        sd_values=sqrt.(1 ./(terminal_nodes_size .*tau .+tau_mu))
        mu_array=randn(n_terminal_node) .*sd_values .+mean_values

        #print(mu_array,"\n")

        for i in 1:length(terminal_nodes_index)
            new_tree.nodes[terminal_nodes_index[i]].mu=mu_array[i]
        end

        return new_tree

    end
end

# Updating the tau value

function update_tau(S,nu,lambda,n)

    # Updating tau value
    tau= rand(Gamma((n+nu)/2,  # Shape parameter
                    ((S+nu*lambda)/2)^(-1) ), # Rate Parameter (th's why it is the inverse)
         1)

    return tau
end

# Update tree prior

function tree_prior(tree_object::tree,alpha::Float64,beta::Float64)

    # Initilaizing the log_tree_prior
    log_tree_prior=0

    # Probability over the root node
    if(length(tree_object.nodes)==0)
            log_tree_prior+=log(1-alpha)

        # Prob. of being not-terminal
    else
            log_tree_prior+=log(alpha)

    end

    # Considering the other nodes

    for i in tree_object.nodes

        # Prob. of being terminal
        if(i.terminal)
            log_tree_prior+=log(1-alpha*(1+i.depth)^(-beta))

        # Prob. of being not-terminal
        else
            log_tree_prior+=log(alpha)-beta*log(1+i.depth)

        end
    end

    return log_tree_prior
end

# Get prediction
function get_prediction(multiple_trees_object::Array,x,single_tree=true)

    # Verify if it is a single tree
    if(length(multiple_trees_object)==1)
        single_tree=true
    end

    # Count the numbr of trees
    n_trees=length(multiple_trees_object)

    # Create an empty vector of the predictions
    predictions=repeat([0.0],size(x)[1])

    # Select just one tree
    if(single_tree)

            if(length(multiple_trees_object[1].nodes)==0)
                predictions.=multiple_trees_object[1].root.mu
            else
                # Getting the terminal nodes index
                for node in multiple_trees_object[1].nodes
                    if(node.terminal)
                            predictions[node.observations_index].=node.mu
                    end
                end
            end
    else
        # Recursion over the other trees
        frac_trees=deepcopy(multiple_trees_object)
        frac_trees=deleteat!(frac_trees,1)
        predictions=get_prediction([multiple_trees_object[1]],x,true)+ #Getting the values of the tree which was blanked
                    get_prediction(frac_trees,x, length(multiple_trees_object)==1)

    end

    return predictions
end

# Main juBART function

function juBart(x::Array, # Matrix of observations
               y::Array{Float64,1}, # Target variable
               number_trees::Int64, # Number of trees
               node_min_size::Int64=5, # Minimum size of observations in each terminal nodes
               alpha::Float64=0.95, # Alpha paramteter from tree prior
               beta::Float64=2.0, # Beta parameter from treee prior
               tau_mu::Float64=1.0, # Tau from mu prior
               nu::Float64=3.0, # Parameter from tau prior
               tau::Float64=0.1, # Initial value to tau
               lambda::Float64=0.1, # Lambda value
               tau_zero::Float64=0.1, # Tau initial value

               n_iter::Int64=1000,
               burn_in::Int64=250,
               thin::Int64=1)

    # Scaling hte y
    y_mean=mean(y)
    y_sd=std(y)
    y_scale=(y.-y_mean) ./y_sd

    sigma=0.0
    # Initializing the algorithm
    n=length(y)
    current_trees=Array{Any}(nothing,number_trees)

    # Think in a error message from number_tree=0

    # Storage containers
    store_size=Int((n_iter-burn_in)/thin)
    tree_store=Array{Any}(nothing,store_size)
    sigma_store=Array{Union{Nothing,Float64},1}(nothing,store_size)
    y_hat_store=Array{Union{Nothing,Float64}}(nothing,store_size,n)
    log_like_store=Array{Union{Nothing,Float64},1}(nothing,store_size)
    full_cond_store=Array{Union{Nothing,Float64},2}(nothing,store_size,number_trees)

    # Creating the stumps
    current_trees=fill(tree(root(x,0,0),tree_node[]),number_trees )

    # Getting the predictions
    predictions=get_prediction(current_trees,x,number_trees==1)
    #print(predictions,"\n")

    # Iterating over MH runs
    @showprogress 1 "MH Runs..." for i in 1:n_iter

        # Creating the storage
        if((i > burn_in) && (i%thin==0))
            curr=convert(Int,round((i-burn_in)/thin,digits=1))
            tree_store[curr]=deepcopy(current_trees)
            sigma_store[curr]=sigma
            y_hat_store[curr,:]=predictions
            #log_like_store[curr]=log_like
        end

        # Iterating over the trees
        for j in 1:number_trees

            # Calculating the residuals for each
            if(number_trees>1)
                partial_trees=deepcopy(current_trees)
                partial_trees=deleteat!(partial_trees,j)

                current_partial_residuals=y_scale-
                                          get_prediction(partial_trees,x, number_trees==1)
            else
                current_partial_residuals=y_scale
            end

            move=rand(["grow","prune","change","swap"],1)[1]
            if i < max(floor(0.1*burn_in), 10) move="grow" end

            #print(move,"\n")
            # Getting the new tree
            new_trees=deepcopy(current_trees)
            new_trees[j]=move_tree(new_trees[j],move,node_min_size)[1]

            #print(new_trees[j].nodes,"\n")
            #print(current_trees[j].root,"\n")

            #print(current_trees,"\n")
            #print(new_trees,"\n")

            #print(move_tree(new_trees[j],move,node_min_size)[1],"\n")

            # Calculating the likelihood of each set
            l_new=tree_full_conditional(new_trees[j],
                                                current_partial_residuals,tau,tau_mu)+
                      tree_prior(new_trees[j], alpha, beta)

            l_old=tree_full_conditional(current_trees[j],
                                            current_partial_residuals,tau,tau_mu)+
                  tree_prior(current_trees[j], alpha, beta)

                  #print("L_old: ",l_old,"\n")
                  #print("L_new: ",l_new,"\n")
            #Probability of accept the new proposed tree
            acceptance=exp(l_new-l_old)
            #print(acceptance,"\n")
            #If Storage or not based on thin and burn parameters

            if((i > burn_in) && ((i %thin) == 0) )
                    full_cond_store[curr, j] = l_old
            end

            if(rand(1)[1]<acceptance)
                    #print("ACCEPT! \n")
                    # Make changes if accept
                    #print(new_trees,"\n")
                    #print(current_trees,"\n")
                    current_trees = deepcopy(new_trees)
                    #print(current_trees,"\n")
                    #print(new_trees,"\n")
            end # End of accept if statement

            #To update the mu values
            current_trees[j]=update_mu(current_trees[j],
                                           current_partial_residuals,
                                           tau,tau_mu)

        end #End of iterations over the trees

     # Get the predicitions from updated trees
     predictions=get_prediction(current_trees,x,number_trees==1)
     S=sum((y_scale-predictions).^2)

    # Update tau
    tau=update_tau(S,nu,lambda,n)[1]

    sigma=1/sqrt(tau)
    #print(sigma)

    # NEED TO THINK HOW TO CALCULATE IT USING JULIA
    #log_like=loglikelihood(Normal(predictions,sigma),y_scale)

    end # End the MH iterations

    return Dict("trees"=>tree_store,
            "sigma"=>sigma_store,
            "predictions"=>y_hat_store,
            "full_conditionals"=>full_cond_store,
            "y"=>y,"X"=>x,"n_iter"=>n_iter,"burn_in"=>burn_in,
            "thin"=>thin,"store_size"=>store_size,"number_trees"=>number_trees)
end

# Function to update a tree with a new X
function predict_tree_model(tree_object::tree,x_new::Array{Float64,2})

    # Creating the new tree
    new_tree=deepcopy(tree_object)

    # Updating the root node
    new_tree.root.x=x_new

    # Getting the index of nodes to be updated
    nodes_index=collect(1:length(new_tree.nodes))

    # Getting the new tree
    new_tree=update_tree_function(new_tree,nodes_index,0)[1]

    return new_tree
end

function predict_juBart(jubart_model::Dict,x_predict::Array{Float64,2},type::String="mean")

    new_jubart_model=deepcopy(jubart_model)
    n_iter=size(new_jubart_model["predictions"])[1]
    n_xpredict=size(x_predict)[1]
    y_hat_matrix=Array{Union{Nothing,Float64}}(nothing,n_iter,n_xpredict)

    @showprogress 1 "Predicting..." for i in 1:n_iter
    #for i in 1:n_iter

        for k in 1:jubart_model["number_trees"]
             new_jubart_model["trees"][i][k]=predict_tree_model(new_jubart_model["trees"][i][k],x_predict)
        end

        y_hat_matrix[i,:]=get_prediction(new_jubart_model["trees"][i],x_predict,jubart_model["number_trees"]==1)

    end

    if(type=="mean")
        return mean(y_hat_matrix,dims=1)
    elseif(type=="median")
        return median(y_hat_matrix,dims=1)
    else
        return y_hat_matrix
    end
end

function sim_friedman_simple(n,p=0,scale_err=1)
    # Simulate some data using a simple version of friedman
    # y= 10sin(πx1*x2)+20(x3=0.5)2+10*4+ 5x5 +ε

    X=rand(n,5+p)
    pars=[10,20,10,5]

    mean=pars[1] .* sin.(pi .*X[:,1] .*X[:,2]) .+ pars[2] .*(X[:,3] .-0.5) .^2 .+
         pars[3] .*X[:,4] .+ pars[4] .*X[:,5]

    y=randn(n).*scale_err .+mean

    return Dict("x"=>X,"y"=>y,"true_mean"=>mean,"true_scale"=>scale_err)
end


text:
module TestStencil
using AnisotropicDiffusion
using AnisotropicDiffusion: laplace, dot_grad
using Test

@testset "laplace" begin
    x = randn(3,3,3)
    index = CartesianIndex((2,2,2))
    @test laplace(x, index) ≈ x[1,2,2] + x[2,1,2] + x[2,2,1] + x[3,2,2] + x[2,3,2] + x[2,2,3] - 6x[2,2,2]

end

@testset "dot_grad" begin
    x = randn(3,3,3)
    y = randn(3,3,3)
    index = CartesianIndex((2,2,2))
    @test dot_grad(x,y,index) ≈
        (x[3,2,2] - x[2,2,2])*(y[3,2,2] - y[2,2,2]) +
        (x[2,3,2] - x[2,2,2])*(y[2,3,2] - y[2,2,2]) +
        (x[2,2,3] - x[2,2,2])*(y[2,2,3] - y[2,2,2])
end

end#module


text:
using Clang
using Clang.LibClang.Clang_jll
pat = raw"C:\Projects\libiio"

# LIBCLANG_HEADERS are those headers to be wrapped.
const LIBCLANG_INCLUDE = joinpath(dirname(Clang_jll.libclang_path), "..", "include", "clang-c") |> normpath
const LIBCLANG_HEADERS = [joinpath(pat, header) for header in readdir(pat) if endswith(header, ".h")]

wc = init(; headers = LIBCLANG_HEADERS,
            output_file = joinpath(@__DIR__, "libiio_api.jl"),
            common_file = joinpath(@__DIR__, "libiio_common.jl"),
            clang_includes = vcat(LIBCLANG_INCLUDE, CLANG_INCLUDE),
            clang_args = ["-I", joinpath(LIBCLANG_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"libiio",
            clang_diagnostics = true,
            )

run(wc)


text:
abstract AbstractRemoteSegments

remotesegments(p::AbstractPolicy,m::AbstractModel,s::AbstractSampler,nproposals::Int) = _remotesegments(p,m,s,nproposals)

numsegments(s::AbstractRemoteSegments) = s.numsegments
numproposalspersegment(s::AbstractRemoteSegments) = s.numproposalspersegment
numtotalproposals(s::AbstractRemoteSegments) = s.numsegments*s.numproposalspersegment

@inline _numjobsegments(::Type{Val{:procs}}) = nprocs()
@inline _numjobsegments(::Type{Val{:workers}}) = nworkers()
@inline _numjobsegments(::Type{Val{:test}}) = 3

@inline _processnumbers(::Type{Val{:procs}}) = procs()
@inline _processnumbers(::Type{Val{:workers}}) = workers()
@inline _processnumbers(::Type{Val{:test}}) = workers()

@inline _numproposalspersegment(nproposals::Int,njobsegments::Int) = ceil(Int,nproposals/njobsegments)


text:
using TFIMSampler

L = parse(Int, ARGS[1])
h = parse(Float64, ARGS[2])
N = parse(Int, ARGS[3])

println("L is $L")
println("h is $h")
println("N is $N")

# pre-compile once

# make directory if none exists
dir = "../april_17_data/"
mkpath(dir)
file = dir * "L-$L-h-$h"

@time configs = TFIMSampler.sample(L=2, h=1.0, N=4, file=file)

@time configs = TFIMSampler.sample(L=L, h=h, N=N, file=file)


text:
# integration with LinearAlgebra stdlib

## transpose and adjoint

function transpose_f!(f, At::AbstractGPUArray{T, 2}, A::AbstractGPUArray{T, 2}) where T
    gpu_call(At, A) do ctx, At, A
        idx = @cartesianidx A
        @inbounds At[idx[2], idx[1]] = f(A[idx[1], idx[2]])
        return
    end
    At
end

LinearAlgebra.transpose!(At::AbstractGPUArray, A::AbstractGPUArray) = transpose_f!(transpose, At, A)
LinearAlgebra.adjoint!(At::AbstractGPUArray, A::AbstractGPUArray) = transpose_f!(adjoint, At, A)

function Base.copyto!(A::AbstractGPUArray{T,N}, B::Adjoint{T, <: AbstractGPUArray{T,N}}) where {T,N}
    adjoint!(A, B.parent)
end

function Base.copyto!(A::AbstractGPUArray{T,N}, B::Transpose{T, <: AbstractGPUArray{T,N}}) where {T,N}
    transpose!(A, B.parent)
end

function Base.copyto!(A::Array{T,N}, B::Adjoint{T, <:AbstractGPUArray{T,N}}) where {T,N}
    copyto!(A, Adjoint(Array(parent(B))))
end
function Base.copyto!(A::Array{T,N}, B::Transpose{T, <:AbstractGPUArray{T,N}}) where {T,N}
    copyto!(A, Transpose(Array(parent(B))))
end

## copy upper triangle to lower and vice versa

function LinearAlgebra.copytri!(A::AbstractGPUMatrix{T}, uplo::AbstractChar) where T
  n = LinearAlgebra.checksquare(A)
  if uplo == 'U'
      gpu_call(A) do ctx, _A
        I = @cartesianidx _A
        i, j = Tuple(I)
        if j > i
          _A[j,i] = _A[i,j]
        end
        return
      end
  elseif uplo == 'L'
      gpu_call(A) do ctx, _A
        I = @cartesianidx _A
        i, j = Tuple(I)
        if j > i
          _A[i,j] = _A[j,i]
        end
        return
      end
  else
      throw(ArgumentError("uplo argument must be 'U' (upper) or 'L' (lower), got $uplo"))
  end
  A
end

## triangular

# mixed CPU/GPU: B -> A
Base.copyto!(A::Array{T,N}, B::UpperTriangular{T, <:AbstractGPUArray{T,N}}) where {T,N} = copyto!(A, UpperTriangular(Array(parent(B))))
Base.copyto!(A::Array{T,N}, B::LowerTriangular{T, <:AbstractGPUArray{T,N}}) where {T,N} = copyto!(A, LowerTriangular(Array(parent(B))))

# GPU/GPU: B -> A
Base.copyto!(A::AbstractGPUArray{T,N}, B::UpperTriangular{T, <:AbstractGPUArray{T,N}}) where {T,N} = LinearAlgebra.triu!(copyto!(A, parent(B)))
Base.copyto!(A::AbstractGPUArray{T,N}, B::LowerTriangular{T, <:AbstractGPUArray{T,N}}) where {T,N} = LinearAlgebra.tril!(copyto!(A, parent(B)))
for T in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular)
    @eval Base.copyto!(A::$T{T, <:AbstractGPUArray{T,N}}, B::$T{T, <:AbstractGPUArray{T,N}}) where {T,N} = $T(copyto!(parent(A), parent(B)))
end

function LinearAlgebra.tril!(A::AbstractGPUMatrix{T}, d::Integer = 0) where T
  gpu_call(A, d; name="tril!") do ctx, _A, _d
    I = @cartesianidx _A
    i, j = Tuple(I)
    if i < j - _d
      _A[i, j] = 0
    end
    return
  end
  return A
end

function LinearAlgebra.triu!(A::AbstractGPUMatrix{T}, d::Integer = 0) where T
  gpu_call(A, d; name="triu!") do ctx, _A, _d
    I = @cartesianidx _A
    i, j = Tuple(I)
    if j < i + _d
      _A[i, j] = 0
    end
    return
  end
  return A
end

## diagonal

Base.copy(D::Diagonal{T, <:AbstractGPUArray{T, N}}) where {T, N} = Diagonal(copy(D.diag))

## matrix multiplication

function generic_matmatmul!(C::AbstractArray{R}, A::AbstractArray{T}, B::AbstractArray{S}, a::Number, b::Number) where {T,S,R}
    if size(A,2) != size(B,1)
        throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $(size(B))"))
    end
    if size(C,1) != size(A,1) || size(C,2) != size(B,2)
        throw(DimensionMismatch("result C has dimensions $(size(C)), needs $((size(A,1),size(B,2)))"))
    end
    if isempty(A) || isempty(B)
        return fill!(C, zero(R))
    end

    gpu_call(C, A, B; name="matmatmul!") do ctx, C, A, B
        idx = @linearidx C
        i, j = @inbounds Tuple(CartesianIndices(C)[idx])..., 1

        @inbounds if i <= size(A,1) && j <= size(B,2)
            z2 = zero(A[i, 1]*B[1, j] + A[i, 1]*B[1, j])
            Ctmp = convert(promote_type(R, typeof(z2)), z2)
            for k in 1:size(A,2)
                Ctmp += A[i, k]*B[k, j]
            end
            C[i,j] = Ctmp*a + C[i,j]*b
        end

        return
    end

    C
end

LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::AbstractGPUVecOrMat, B::AbstractGPUVecOrMat, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::AbstractGPUVecOrMat, B::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::AbstractGPUVecOrMat, B::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, B::AbstractGPUVecOrMat, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, B::AbstractGPUVecOrMat, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, B::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, B::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, B::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, B::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, a::Number, b::Number) = generic_matmatmul!(C, A, B, a, b)

# specificity hacks
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::AbstractGPUVecOrMat, B::AbstractGPUVecOrMat, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::AbstractGPUVecOrMat, B::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::AbstractGPUVecOrMat, B::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, B::AbstractGPUVecOrMat, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, B::AbstractGPUVecOrMat, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, B::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, B::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, B::LinearAlgebra.Adjoint{<:Any, <:AbstractGPUVecOrMat}, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)
LinearAlgebra.mul!(C::AbstractGPUVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, B::LinearAlgebra.Transpose{<:Any, <:AbstractGPUVecOrMat}, a::Real, b::Real) = generic_matmatmul!(C, A, B, a, b)

function generic_rmul!(X::AbstractArray, s::Number)
    gpu_call(X, s; name="rmul!") do ctx, X, s
        i = @linearidx X
        @inbounds X[i] *= s
        return
    end
    return X
end

LinearAlgebra.rmul!(A::AbstractGPUArray, b::Number) = generic_rmul!(A, b)

function generic_lmul!(s::Number, X::AbstractArray)
    gpu_call(X, s; name="lmul!") do ctx, X, s
        i = @linearidx X
        @inbounds X[i] = s*X[i]
        return
    end
    return X
end

LinearAlgebra.lmul!(a::Number, B::AbstractGPUArray) = generic_lmul!(a, B)

## permutedims

function LinearAlgebra.permutedims!(dest::AbstractGPUArray, src::AbstractGPUArray,
                                    perm::NTuple)
    Base.checkdims_perm(dest, src, perm)
    function permutedims_kernel(ctx, dest, src, ::Val{perm}) where {perm}
        I = @cartesianidx src
        @inbounds begin
            J = CartesianIndex(map(i->I[i], perm))
            dest[J] = src[I]
        end
        return
    end
    gpu_call(permutedims_kernel, dest, src, Val(perm))
    return dest
end

# TODO: implementation without the memory copy
LinearAlgebra.permutedims!(dest::AbstractGPUArray, src::AbstractGPUArray, perm) =
    permutedims!(dest, src, Tuple(perm))

## norm

function LinearAlgebra.norm(v::AbstractGPUArray{T}, p::Real=2) where {T}
    norm_x = if p == Inf
        maximum(abs.(v))
    elseif p == -Inf
        minimum(abs.(v))
    else
        mapreduce(x->abs(x)^p, +, v; init=float(zero(T)))^(1/p)
    end
    return real(norm_x)
end

## symmetric

# prevent scalar indexing (upstream? this version is slower than a simple loop)
function Base.similar(A::Hermitian{<:Any,<:AbstractGPUArray}, ::Type{T}) where T
    B = similar(parent(A), T)
    fill!(view(B, diagind(B)), 0)
    return Hermitian(B, ifelse(A.uplo == 'U', :U, :L))
end


text:
export plot_equidistant_2d_mesh

# Shifts a regularly spaced list of coordinates
function center_to_corner_coords(xs)
    @assert length(xs) > 0 "Need at least one point to plot per axis"
    if length(xs) == 1
        # Plot a single value as
        [-0.5, 0.5] .+ xs[1]
    else
        step = xs[2] - xs[1]
        [xs; [xs[end] + step]] .- step / 2
    end
end

"""
Generates a 2D plot where values are represented by uniformly sized coloured
rectangles.

Parameters
----------
xs
    A one-dimensional array of the data point x coordinates (centered). Assumed
    to be equally spaced.
ys
    A one-dimensional array of the data point y coordinates (centered). Assumed
    to be equally spaced.
zs
    A two-dimensional array of the data values at the corresponding coordinates.
    Outer dimension corresponds to the x axis.
ax
    The Axes object to use. Defaults to gca().

Returns
-------
A tuple with the Axes and plot objects.
"""
function plot_equidistant_2d_mesh(xs::Vector, ys::Vector, zs::Matrix, ax=nothing; kwargs...)
    if ax === nothing
        ax = gca()
    end

    corner_xs = center_to_corner_coords(xs)
    corner_ys = center_to_corner_coords(ys)
    img = ax.pcolormesh(corner_xs, corner_ys, zs'; kwargs...)
    ax.set_xlim(corner_xs[1], corner_xs[end])
    ax.set_ylim(corner_ys[1], corner_ys[end])

    ax, img
end


text:
using HDF5

data_dir = "../data/modelnet40_ply_hdf5_2048"

cd(@__DIR__)
test_files = readlines("../data/modelnet40_ply_hdf5_2048/test_files.txt")
train_files = readlines("../data/modelnet40_ply_hdf5_2048/train_files.txt")


text:
#Written By Johnathan Bizzano

function dense_kernel(inputs, constant_row, activ, idx)
    sum = 0.0
    for i in 1:row_len
        sum += inputs[i] * constant_row[idx + i]
    end
    activ(sum)
end

function Dense(size, initializer; activation=Activation.None())
    return LayerGenerator(
        function(pb, sett, input_shape)
                indexes = prealloc_readonly!(pb, [(i - 1) * size for i in 1:size])
                activ = activation(sett)
                device = sett.device

                LayerDesign(Sig(:DenseLayer), input_shape, (size), (size, size), initializer,
                        LayerInitializer(
                            function (inputs, outputs, constants)
                                indexes = indexes()
                                inputs = Ref(inputs)
                                constants = Ref(constants)
                                return () -> broadcast!(dense_kernel, outputs, inputs, constants, activ, size, indexes)
                            end))
    end)
end

export DenseLayer
DenseLayer(size, initializer; activation=Activation.None()) = Dense(size, initializer; activation=activation)


text:
export SeparableFcn, AddedSeparableFcn, ui,
 uk, σscale, D!, D, H!, H

abstract type SeparableFcn end

## Need tools to identify

struct AddedSeparableFcn{T <: Tuple} <: SeparableFcn
    list::T
end

function Base.show(io::IO, Σg::AddedSeparableFcn)
    println(io, "AddedSeparableFcn:")
    for g in Σg.list
        println(io, "  $g")
    end
end

# """
#     g₁::SeparableFcn + g₂::SeparableFcn
#
# Add the separable functions so that `(g₁ + g₂)(z) = g₁(z) + g₂(z)`.
# """

+(g::SeparableFcn, Σg::AddedSeparableFcn) = AddedSeparableFcn((g, Σg.list...))
+(Σg::AddedSeparableFcn, g::SeparableFcn) = AddedSeparableFcn((Σg.list..., g))
function +(Σg₁::AddedSeparableFcn, Σg₂::AddedSeparableFcn)
    AddedSeparableFcn((Σg₁..., Σg₂...))
end

+(g::SeparableFcn...) = AddedSeparableFcn(g)

(Σg::AddedSeparableFcn)(z) = mapreduce(gi->gi(z),+, Σg.list)

# (+)(f1::T, f2::U) where {T<:SeparableFcn, U<:SeparableFcn}  =

### Generate parameterization for the i-th component of a lower triangular map U
#
#
#  U: z =(z1, z2,...,zn) → [U1(z1);
#                          U2(z1, z2);
#                            ...
#                          Un(z1, z2, ..., zn)]
#
#

# We use separable maps RadialMapComponent(z1,...,zk) = Σi=1,k-1 ui(z1,..,zi) + uk(z1,..,zk)
# uk has a different structure from ui to ensure the monotonocity condition
# of the triangular map
struct ui <:SeparableFcn
        p::Int64
        ξi::Array{Float64,1}
        σi::Array{Float64,1}
        coeffi::Array{Float64,1}
        # f::Array{ParamFcn,1}

        # Inner constructor
        function ui(p::Int64, ξi::Array{Float64, 1}, σi::Array{Float64, 1}, coeffi::Array{Float64, 1})

                if p==-1
                        # p==-1 , ui(z) = 0.0
                        @assert size(ξi,1)==0 "Size of ξi does not match the number of RBFs p"
                        @assert size(σi,1)==0 "Size of σi does not match the number of RBFs p"
                        @assert size(coeffi,1)==0 "Size of coeffi does not match the number of RBFs p"
                else
                        # p==0 ui(z) = az, p>0 ui(z) = a1z + p rbfs
                        @assert size(ξi,1)==p "Size of ξi does not match the number of RBFs p"
                        @assert size(σi,1)==p "Size of σi does not match the number of RBFs p"
                        @assert size(coeffi,1)==p+1 "Size of coeffi does not match the number of RBFs p"
                end
                new(p, ξi, σi, coeffi)
        end
end

function ui(p::Int64)
        if p==-1
                return ui(p, Float64[], Float64[], Float64[])
        else
                return ui(p, zeros(p), ones(p), zeros(p+1))
        end
end

function (u::ui)(z::Real)
        p = u.p
if p==-1
        return 0.0
else
        out = u.coeffi[1]*z
        if p>0
                @inbounds for j=1:p
                out +=rbf(u.ξi[j], u.σi[j])(z)*u.coeffi[j+1]
                end
        end
        return out
end
end

# Compute derivative of a ui function
function D!(u::ui, z::Real)
        p = u.p
if p==-1
        return 0.0
else
        out = u.coeffi[1]
        if p>0
                for j=1:p
                @inbounds out +=rbf′(u.ξi[j], u.σi[j])(z)*u.coeffi[j+1]
                end
        end
        return out
end
end

D(u::ui) =  z-> D!(u, z)

#### Structure for the final component of RadialMapComponent

# uk(z) = c + Σj=0, p+1 ukj ψⱼ(z), p.36 Couplings for nonlinear ensemble filtering

struct uk <:SeparableFcn
        p::Int64
        ξk::Array{Float64,1}
        σk::Array{Float64,1}
        coeffk::Array{Float64,1}

        function uk(p::Int64, ξk::Array{Float64, 1}, σk::Array{Float64, 1}, coeffk::Array{Float64, 1})
                if p==-1
                @assert size(ξk,1)==0 "Size of ξk doesn't match the number of RBFs p"
                @assert size(σk,1)==0 "Size of σk doesn't match the number of RBFs p"
                @assert size(coeffk,1)==0 "Size of coeffk doesn't match the number of RBFs p"
                elseif p==0
                @assert size(ξk,1)==0 "Size of ξk doesn't match the number of RBFs p"
                @assert size(σk,1)==0 "Size of σk doesn't match the number of RBFs p"
                @assert size(coeffk,1)==2 "Size of coeffk doesn't match the number of RBFs p"
                else
                @assert size(ξk,1)==p+2 "Size of ξk doesn't match the number of RBFs p"
                @assert size(σk,1)==p+2 "Size of σk doesn't match the number of RBFs p"
                @assert size(coeffk,1)==p+3 "Size of coeffk doesn't match the number of RBFs p"
                end
                return new(p, ξk, σk, coeffk)
        end
end

# uk is an affine function if p=0

function uk(p::Int64)
        if p==-1
                # uk(z) = z
                return uk(p, Float64[], Float64[], Float64[])
        elseif p==0
                # uk(z) = a[1] + a[2]*z
                return uk(p, Float64[], Float64[], zeros(2))
        else
                return uk(p, zeros(p+2), ones(p+2), zeros(p+3))
        end
end

function (u::uk)(z::T) where {T<:Real}
        p = u.p
        if p==-1
                return z
        elseif p==0
                out = u.coeffk[1] + u.coeffk[2]*z
                return out
        else
                out = u.coeffk[1]
                out +=ψ₀(u.ξk[1], u.σk[1])(z)*u.coeffk[2]
                out +=ψpp1(u.ξk[p+2], u.σk[p+2])(z)*u.coeffk[p+3]

                @inbounds for j=2:p+1
                out +=ψj(u.ξk[j], u.σk[j])(z)*u.coeffk[j+1]
                end
                return out
        end
end

# Compute derivative of a uk function
function D!(u::uk, z::Real)
        p = u.p

        if p==-1
                return 1.0
        elseif p==0
        # Linear function
                return u.coeffk[2]
        else
                out = ψ₀′(u.ξk[1], u.σk[1])(z)*u.coeffk[2]
                out += ψpp1′(u.ξk[p+2], u.σk[p+2])(z)*u.coeffk[p+3]

                @inbounds for j=2:p+1
                out +=rbf(u.ξk[j], u.σk[j])(z)*u.coeffk[j+1]
                end
                return out
        end
end

D(u::uk) =  z-> D!(u, z)

# Compute hessian of a uk function
function H!(u::uk, z::Real)
        p = u.p
        if p<1
                return 0.0
        else
                # ψ₀'' = -rbf
                out =rbf(u.ξk[1], u.σk[1])(z)*(-u.coeffk[2])
                # ψpp1′′
                out +=rbf(u.ξk[p+2], u.σk[p+2])(z)*u.coeffk[p+3]

                @inbounds for j=2:p+1
                out +=rbf′(u.ξk[j], u.σk[j])(z)*u.coeffk[j+1]
                end
                return out
        end
end

H(u::uk) =  z-> H!(u, z)

function σscale(ξ::Array{Float64,1}, γ::Float64)
        p = size(ξ,1)
        σ = zero(ξ)
        # Special treatment for the edges
        # ξ[0]=ξ[1] & ξ[p+1]=ξ[p]
        σ[1] = ξ[2]-ξ[1]
        σ[end] = ξ[end]-ξ[end-1]
        for j=2:p-1
        @inbounds σ[j] = 0.5*(ξ[j+1]-ξ[j-1])
        end

        rmul!(σ, γ)
        return σ
end


text:
# Autogenerated wrapper script for Htop_jll for armv6l-linux-musleabihf
export htop

using Ncurses_jll
JLLWrappers.@generate_wrapper_header("Htop")
JLLWrappers.@declare_executable_product(htop)
function __init__()
    JLLWrappers.@generate_init_header(Ncurses_jll)
    JLLWrappers.@init_executable_product(
        htop,
        "bin/htop",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()