text:
# Feb 2016, orthopteroid@gmail.com, MIT license.
# for julia 0.5.0-dev+2944 (e01b357)
# kudos to https://github.com/rennis250/SDL.jl

macro checked_lib(libname, path)
    ( Base.Libdl.dlopen_e(path) == C_NULL ) && error("Unable to load $libname ($path).")
    quote const $(esc(libname)) = $path end
end
@checked_lib libSDL "/usr/lib/x86_64-linux-gnu/libSDL-1.2.so.0"

type SDL_AudioSpec
        freq::Int32
        format::UInt16
        channels::UInt8
        silence::UInt8
        samples::UInt16
        padding::UInt16
        size::UInt32
		callback::Ptr{Void}
        userdata::Ptr{Void}
end

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

# represents current phase angle and direction of sweep
type Oscillator
	theta::Cfloat
	step::Cfloat
end

Oscillator(freq::Real, sampledelta_1khz::Real) = Oscillator( Cfloat(0), Cfloat(freq * sampledelta_1khz) )

# step the oscillator to the next phase angle for the programmed frequency
@inline function step!(O::Oscillator)
	O.theta += O.step
	if O.theta < Cfloat(-1.57079) # pi / 2
		O.theta -= (O.theta - Cfloat(-1.57079))
		O.step = -O.step
	elseif O.theta > Cfloat(+1.57079)
		O.theta -= (O.theta - Cfloat(+1.57079))
		O.step = -O.step
	end
	return O
end

# given the oscillator's angle, return sin(angle)
@inline function signal(O::Oscillator)
	# only valid between -pi/2 and +pi/2
	const Cpi = Cfloat(3.14159)
	y = ((Cfloat(4.) / Cpi) * O.theta + (Cfloat(-4.) / (Cpi * Cpi)) * O.theta * abs( O.theta ))
	z = Cfloat(0.225) * (y * abs( y ) - y) + y
	return Cshort( round( z * 32765 ) )
end

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

const SAMPLES = 4096 * 1 # bigger buffer = less skipping
const SAMPLERATE = 22050
const SAMPLEDELTA_1KHZ = 2 * pi / SAMPLERATE

# lame attempt at double buffering
const buffer0 = zeros(Cshort, SAMPLES)
const buffer1 = zeros(Cshort, SAMPLES)
const bufbool = Array(Bool)

# this darn thing skips - I think julia is taking too long in this callpath
function mixerCallback(userdata::Ptr{Void}, stream::Ptr{Cshort}, len::Int32)
	rc = ccall((:SDL_MixAudio, libSDL), Int32, (Ptr{Cshort}, Ref{Cshort}, UInt32, Int32), stream, bufbool[] ? buffer1 : buffer0, len, 0x80) # SDL_MIX_MAXVOLUME
	return nothing
end
const mixerCallback_c = cfunction(mixerCallback, Void, (Ptr{Void}, Ptr{Cshort}, Int32))

asDesired = SDL_AudioSpec( SAMPLERATE, 0x9010, 1, 0, SAMPLES, 0, 0, mixerCallback_c, 0 ) # AUDIO_S16MSB
asCurrent = SDL_AudioSpec( 0, 0, 0, 0, 0, 0, 0, mixerCallback_c, 0 )

function main()
	global bufbool # hack: redeclare as a global
	bufbool[] = false

	rc = ccall((:SDL_Init, libSDL), Int32, (UInt32, ), 0x0001 | 0x0010 ) # SDL_INIT_TIMER | SDL_INIT_AUDIO
	rc = ccall((:SDL_OpenAudio, libSDL), Int32, (Ref{SDL_AudioSpec}, Ref{SDL_AudioSpec}), asDesired, asCurrent)
	rc = ccall((:SDL_PauseAudio, libSDL), Void, (UInt32, ), 0x00 ) # arg is pause_on

	osc = Oscillator( 440, SAMPLEDELTA_1KHZ )

	lasttick = ccall((:SDL_GetTicks, libSDL), UInt32, ()) + 1 * 500
	while ccall((:SDL_GetTicks, libSDL), UInt32, ()) <= lasttick

		# fill the buffer that is not playing
		if bufbool[] == false
			for i = 1:SAMPLES
				buffer1[ i ] = signal( step!( osc ) )
			end
		else
			for i = 1:SAMPLES
				buffer0[ i ] = signal( step!( osc ) )
			end
		end

		# swap playback and fill buffers
		rc = ccall((:SDL_LockAudio, libSDL), Void, ())
		bufbool[] = !bufbool[]
		rc = ccall((:SDL_UnlockAudio, libSDL), Void, ())

	end

	rc = ccall((:SDL_Quit, libSDL), Void, ())
	return 0
end

max_stars_repo_name: orthopteroid/JuliaSDLAudio

max_issues_repo_name: orthopteroid/JuliaSDLAudio

max_forks_repo_name: orthopteroid/JuliaSDLAudio


text:
# fisher information
export fisher_det

# TODO: mega wrong, don't pass in the uav's sensor, pass in the filter's
fisher_det(uav::Vehicle, df::DF) = fisher_det(uav.sensor, df)

fisher_det(df) = fisher_det(df.sensor, df)

function fisher_det(bo::BearingOnly, df::DF)
	fd = zeros(df.n, df.n)
	ns2 = bo.noise_sigma * bo.noise_sigma

    # Loop over possible vehicle states
	for xi = 1:df.n, yi = 1:df.n
        x = (xi - 0.5) * df.cell_size
        y = (yi - 0.5) * df.cell_size

        # Loop over possible target states
        a = bc = d = 0.0
        for txi = 1:df.n, tyi = 1:df.n
            tx = (txi - 0.5) * df.cell_size
            ty = (tyi - 0.5) * df.cell_size
            p = df.b[txi,tyi]

            dx = tx - x
            dy = ty - y

            # ensure dx, dy are non-zero (to avoid singularities)
            dx == 0.0 && (dx = 1e-6)
            dy == 0.0 && (dy = 1e-6)

            r2 = dx*dx + dy*dy
            r4 = r2*r2

            a += p*dy*dy / r4
            bc -= p*dx*dy / r4
            d += p*dx*dx / r4
        end
        a /= ns2
        bc /= ns2
        d /= ns2

        fd[xi, yi] = a*d - bc*bc
	end

	return fd
end

function fisher_det(bo::RangeOnly, df::DF)
	fd = zeros(df.n, df.n)
	ns2 = bo.noise_sigma * bo.noise_sigma
	for xind = 1:df.n
		x = (xind - 0.5) * df.cell_size
		for yind = 1:df.n
			y = (yind - 0.5) * df.cell_size

			# now we must loop over belief
			a = bc = d = 0.0
			for txi = 1:df.n
				tx = (txi - 0.5) * df.cell_size
				for tyi = 1:df.n
					ty = (tyi - 0.5) * df.cell_size
					p = df.b[txi,tyi]
					dx = tx - x
					dy = ty - y
					if dx == 0.0
						dx = 1e-6
					end
					if dy == 0.0
						dy = 1e-6
					end
					r2 = dx*dx + dy*dy

					a += p*dx*dx / r2
					bc += p*dx*dy / r2
					d += p*dy*dy / r2
				end
			end
			a /= ns2
			bc /= ns2
			d /= ns2

			fd[xind, yind] = a*d - bc*bc

		end
	end
	return fd
end

ext: jl


text:
module BDMAnalytic

    export BDM, prob, SSprob

    using Plots, Parameters, LinearAlgebra, Distributions, StatsBase, DoubleFloats

    """
    Define the binary decision model struct:
        * N: the number of agents in the system.
        * pars: the parameters of the system. Note the length of pars specifies the model type, can be [ϵ,μ], [ϵ1,ϵ2,μ] or [ϵ1,ϵ2,μ1,μ2].
    """
    @with_kw struct BDM
        N::Int64 = 100;
        pars::Vector{Float64}; # can be [ϵ,μ], [ϵ1,ϵ2,μ] or [ϵ1,ϵ2,μ1,μ2]
        model_type::String = "ant" # be also be "voter" or "brock"
        A::Matrix{Float64} = make_TRM(pars, N, model_type);
        λ::Array{Complex{Double64}} = GetEigVals(A,pars,N,model_type);
        q_arr::Vector{Vector{Complex{Double64}}} = [GetOrthoQ(pars,N,λ[j],model_type) for j in 1:N+1];
        p_arr::Vector{Vector{Complex{Double64}}} = [GetOrthoP(pars,N,λ[j],model_type) for j in 1:N+1];
        den_prod::Vector{Complex{Double64}} = [prod([λ[i]-λ[j] for j in filter!(e->e≠i,[j for j in 1:N+1])]) for i in 1:N+1];
        As::Vector{Double64} = [a(pars,N,j,model_type) for j in 1:N];
        Bs::Vector{Double64} = [b(pars,N,j,model_type) for j in 0:N-1];
    end

    """
    Define rate function for ants going to the rh food source
    """
    function a(pars::Vector{Float64}, N::Int64, n::Int64, model_type::String)
        if model_type == "ant"
            if length(pars)==2
                ϵ = pars[1];
                μ = pars[2];
                return convert(Double64,(N-(n-1))*ϵ+μ*(n-1)*(N-(n-1))/(N-1))::Double64
            elseif length(pars)==3
                ϵ1 = pars[1];
                μ = pars[3];
                return convert(Double64,(N-(n-1))*ϵ1+μ*(n-1)*(N-(n-1))/(N-1))::Double64
            elseif length(pars)==4
                ϵ1 = pars[1];
                μ1 = pars[3];
                return convert(Double64,(N-(n-1))*ϵ1+μ1*(n-1)*(N-(n-1))/(N-1))::Double64
            else
                error("length of pars must be 2, 3 or 4.")
            end
        elseif model_type == "voter"
            if length(pars)==1
                p₀ = pars[1];
                if p₀ > 1 || p₀ < 0
                    error("p₀ must be between 0 and 1.")
                end
                ϵ = p₀/(2*N);
                μ = (1-p₀)*(N-1)/(2*N^2);
                return convert(Double64,(N-(n-1))*ϵ+μ*(n-1)*(N-(n-1))/(N-1))::Double64
            elseif length(pars)==2
                ϵ = pars[1];
                μ = pars[2];
                return convert(Double64,(N-(n-1))*ϵ+μ*(n-1)*(N-(n-1))*(1+(N-(n-1))/(N-1))/(N-1))::Double64
            else
                error("length of pars must be 1 or 2 for the voter models")
            end
        elseif model_type == "brock_exact" || model_type == "brock"
            if length(pars)==5
                α = pars[1];
                β = pars[2];
                γ = pars[3];
                F = pars[4];
                J = pars[5];
                mn = (2*(n-1)-N)/N;
                return convert(Double64,(N-(n-1))*γ/(1+exp(-β*(F+J*(α+1)*mn))))
            else
                error("length of pars must be 5 for the brock and durlauf model")
            end
        else
            error("chosen model type is not recognised")
        end
    end

    """
    Define rate function for ants going to the lh food source
    """
    function b(pars::Vector{Float64}, N::Int64, n::Int64, model_type::String)
        if model_type == "ant"
            if length(pars)==2
                ϵ = pars[1];
                μ = pars[2];
                return convert(Double64,(n+1)*ϵ+μ*(n+1)*(N-(n+1))/(N-1))::Double64
            elseif length(pars)==3
                ϵ2 = pars[2];
                μ = pars[3];
                return convert(Double64,(n+1)*ϵ2+μ*(n+1)*(N-(n+1))/(N-1))::Double64
            elseif length(pars)==4
                ϵ2 = pars[2];
                μ2 = pars[4];
                return convert(Double64,(n+1)*ϵ2+μ2*(n+1)*(N-(n+1))/(N-1))::Double64
            else
                error("length of pars must be 2, 3 or 4.")
            end
        elseif model_type == "voter"
            if length(pars)==1
                p₀ = pars[1];
                ϵ = p₀/(2*N);
                μ = (1-p₀)*(N-1)/(2*N^2);
                return convert(Double64,(n+1)*ϵ+μ*(n+1)*(N-(n+1))/(N-1))::Double64
            elseif length(pars)==2
                ϵ = pars[1];
                μ = pars[2];
                return convert(Double64,(n+1)*ϵ+μ*(n+1)*(N-(n+1))*(1+(n+1)/(N-1))/(N-1))::Double64
            else
                error("length of pars must be 1 or 2 for the voter models")
            end
        elseif model_type == "brock_exact" || model_type == "brock"
            if length(pars)==5
                α = pars[1];
                β = pars[2];
                γ = pars[3];
                F = pars[4];
                J = pars[5];
                mn = (2*(n+1)-N)/N;
                return convert(Double64,(n+1)*γ/(1+exp(β*(F+J*(α+1)*mn))))
            else
                error("length of pars must be 5 for the brock and durlauf model")
            end
        else
            error("chosen model type is not recognised")
        end
    end

    """
    Make the transition rate matrix
    """
    function make_TRM(pars::Vector{Float64}, N::Int64, model_type::String)
        A = zeros(N+1, N+1);
        for i in 1:size(A)[1]
            for j in 1:size(A)[2]
                if i == 1 && j == 1
                    A[1,1] = -a(pars,N,1, model_type)
                elseif i == j && i>1
                    A[i,i] = -(a(pars,N,i, model_type)+b(pars,N,i-2, model_type))
                elseif i == j+1
                    A[i,j] = a(pars,N,j, model_type)
                elseif i == j-1
                    A[i,j] = b(pars,N,i-1, model_type)
                else
                    continue
                end
            end
        end
        return A::Matrix{Float64}
    end

    """
    Get eigenvalues
    """
    function GetEigVals(A::Matrix{Float64},pars::Vector{Float64},N::Int64,model_type::String)
        if model_type == "ant"
            if length(pars)==2
                ϵ = pars[1];
                μ = pars[2];
                λ = [-(m-1)*(2*ϵ+(m-2)μ/(N-1)) for m in 1:N+1];
                return convert(Array{Complex{Double64}},λ)
            elseif length(pars)==3
                ϵ1 = pars[1];
                ϵ2 = pars[2]
                μ = pars[3];
                λ = [-(m-1)*(ϵ1+ϵ2+(m-2)μ/(N-1)) for m in 1:N+1];
                return convert(Array{Complex{Double64}},λ)
            elseif length(pars)==4
                λ = convert(Array{Complex{Double64}}, reverse(eigvals(A)));
                if λ[1] == λ[2] # if get repeated zero eigenvalues from solver manually separate.
                    λ[1] = 0.0 + (1E-30)im; λ[2] = 0.0 - (1E-30)im;
                end
                return λ::Array{Complex{Double64}}
            else
                error("length of pars must be 2, 3 or 4.")
            end
        elseif model_type == "voter"
            if length(pars)==1
                p₀ = pars[1];
                ϵ = p₀/(2*N);
                μ = (1-p₀)*(N-1)/(2*N^2);
                λ = [-(m-1)*(2*ϵ+(m-2)μ/(N-1)) for m in 1:N+1];
                return convert(Array{Complex{Double64}},λ)
            elseif length(pars)==2
                λ = convert(Array{Complex{Double64}}, reverse(eigvals(A)));
                if λ[1] == λ[2] # if get repeated zero eigenvalues from solver manually separate.
                    λ[1] = 0.0 + (1E-30)im; λ[2] = 0.0 - (1E-30)im;
                end
                return λ::Array{Complex{Double64}}
            else
                error("length of pars must be 1 or 2 for voter models")
            end
        elseif model_type == "brock"
            if length(pars)==5
                α = pars[1];
                β = pars[2];
                γ = pars[3];
                F = pars[4];
                J = pars[5];
                λ = [-(m-1)*(N-β*J*(α+1)*(2+N-m))*γ/N for m in 1:N+1];
                return convert(Array{Complex{Double64}},λ)
            else
                error("length of pars must be 5 for brock model")
            end
        elseif model_type == "brock_exact"
            if length(pars)==5
                λ = convert(Array{Complex{Double64}}, reverse(eigvals(A)));
                if λ[1] == λ[2] # if get repeated zero eigenvalues from solver manually separate.
                    λ[1] = 0.0 + (1E-30)im; λ[2] = 0.0 - (1E-30)im;
                end
                return λ::Array{Complex{Double64}}
            else
                error("length of pars must be 5 for brock model")
            end
        end
    end

    """
    Get the p orthogonal polynomials
    """
    function GetOrthoP(pars::Vector{Float64}, N::Int64, λᵢ::Complex{Double64},model_type::String)
        p = Array{Complex{Double64},1}(undef, N+1);
        p[1] = 1.0; p[2] = λᵢ+a(pars,N,1,model_type);
        for i in 3:N+1
            p[i] = (λᵢ+a(pars,N,i-1,model_type)+b(pars,N,i-3,model_type))*p[i-1] - b(pars,N,i-3,model_type)*a(pars,N,i-2,model_type)*p[i-2]
        end
        return p::Vector{Complex{Double64}}
    end

    """
    Get the q orthogonal polynomials
    """
    function GetOrthoQ(pars::Vector{Float64}, N::Int64, λᵢ::Complex{Double64},model_type::String)
        q = Array{Complex{Double64},1}(undef, N+3);
        q[N+3] = 1.0; q[N+2] = λᵢ + b(pars,N,N-1,model_type);
        for i in reverse([j for j in 3:N+1])
            q[i] = (λᵢ + a(pars,N,i-1,model_type)+b(pars,N,i-3,model_type))*q[i+1] - b(pars,N,i-2,model_type)*a(pars,N,i-1,model_type)*q[i+2]
        end
        return q::Vector{Complex{Double64}}
    end

    """
    Define the sum of the elements from the solution
    """
    function sum_elems(λᵢ::Complex{Double64}, t::Float64, m::Int64, m₀::Int64, pars::Array{Float64,1}, N::Int64, p_arrᵢ::Vector{Complex{Double64}}, q_arrᵢ::Vector{Complex{Double64}}, den_prodᵢ::Complex{Double64})
        return exp(λᵢ*t)*p_arrᵢ[m+1]*q_arrᵢ[m₀+3]/den_prodᵢ::Complex{Double64}
    end

    """
    Define P(x,t|x₀)
    """
    function pm(t::Float64, BD::BDM, m::Int64, m₀::Int64)
        @unpack λ, pars, N, As, Bs, q_arr, p_arr, den_prod = BD
        if m<m₀
            return convert(Double64,N)*prod(Bs[m+1:m₀])*sum([sum_elems(λ[i], t, m, m₀, pars, N, p_arr[i], q_arr[i], den_prod[i]) for i in 1:N+1])::Complex{Double64}
        elseif m==m₀
            return convert(Double64,N)*sum([sum_elems(λ[i], t, m, m, pars, N, p_arr[i], q_arr[i], den_prod[i]) for i in 1:N+1])::Complex{Double64}
        else
            return convert(Double64,N)*prod(As[m₀+1:m])*sum([sum_elems(λ[i], t, m₀, m, pars, N, p_arr[i], q_arr[i], den_prod[i]) for i in 1:N+1])::Complex{Double64}
        end
    end

    """
    Define the probability distribution return function from a initial distribution.
    """
    function prob(BD::BDM, t::Float64, q_init_D::Distribution{Univariate, Discrete})
        @unpack N = BD
        q_init = pdf(q_init_D)
        pmt = Array{Complex{Double64}}(undef,N+1)
        for i in 1:N+1 # loop over the m0's
            pmt[i] = sum([q_init[n+1]*pm(t, BD, i-1, n) for n in 0:N])
        end
        norm = sum(real(pmt))/N
        return (LinRange(0.0,1.0,N+1),real(pmt)/norm)::Tuple{LinRange{Float64}, Vector{Double64}}
    end

    """
    Define the probability distribution return function from a precise value of m₀.
    """
    function prob(BD::BDM, t::Float64, m₀::Int64)
        @unpack N = BD
        norm = sum(real([pm(t, BD, m, m₀) for m in 0:N]))/N
        return (LinRange(0.0,1.0,N+1),real([pm(t, BD, m, m₀) for m in 0:N])/norm)::Tuple{LinRange{Float64}, Vector{Double64}}
    end

    """
    Define the probability distribution return function for N/2 value of m₀.
    """
    function prob(BD::BDM, t::Float64)
        @unpack N = BD
        m₀ = floor(Int64,N/2) # note that we have the n=0 state too.
        norm = sum(real([pm(t, BD, m, m₀) for m in 0:N]))/N
        return (LinRange(0.0,1.0,N+1),real([pm(t, BD, m, m₀) for m in 0:N])/norm)::Tuple{LinRange{Float64}, Vector{Double64}}
    end

    """
    Return the steady state distribution
    """
    function SSprob(BD::BDM)
        @unpack As, Bs, N = BD
        ps = Vector{Double64}(undef, N+1)
        for m in 2:N # use the product rule
            ps[m] = prod(As[1:m-1])*prod(Bs[m:N])
        end
        ps[1] = prod(Bs[1:N]) # do product for the B's
        ps[N+1] = prod(As[1:N]) # do product for the A's
        return (LinRange(0.0,1.0,N+1),N.*ps/sum(ps))::Tuple{LinRange{Float64}, Vector{Double64}}
    end

end # module end

max_stars_repo_name: jamesholehouse/AntModels

max_issues_repo_name: jamesholehouse/AntModels

max_forks_repo_name: jamesholehouse/AntModels


text:
function num_allocs(model; evals=10, samples=2)
    dt = 0.1
    x, u = rand(model)
    n,m = size(model)
    z = KnotPoint(x, u, dt)
    ∇c  = zeros(n,n+m)
    b1 = @benchmark dynamics($model, $x, $u) evals=evals samples=samples
    b2 = @benchmark jacobian!($∇c, $model, $z) evals=evals samples=samples
    b3 = @benchmark discrete_dynamics($RK3, $model, $x, $u, $z.t, $dt) evals=evals samples=samples
    b4 = @benchmark discrete_jacobian!($RK3, $∇c, $model, $z) evals=evals samples=samples
    return b1.allocs + b2.allocs + b3.allocs + b4.allocs
end

# Acrobot
acrobot = RobotZoo.Acrobot()
@test size(acrobot) == (4,1)
@test num_allocs(acrobot) == 0

# Car
car = RobotZoo.DubinsCar()
@test size(car) == (3,2)
@test num_allocs(car) == 0

# Bicycle Car
bicycle = RobotZoo.BicycleModel()
@test size(bicycle) == (4,2)
@test num_allocs(car) == 0

bicycle = RobotZoo.BicycleModel(ref=:rear)
@test size(bicycle) == (4,2)
@test num_allocs(car) == 0

# Planar Rocket
rocket = RobotZoo.PlanarRocket()
@test size(rocket) == (8,2)
@test num_allocs(rocket) == 0

# Planar Quad
quad = RobotZoo.PlanarQuadrotor()
@test size(quad) == (6,2)
@test num_allocs(rocket) == 0

# Cartpole
cartpole = RobotZoo.Cartpole()
@test size(cartpole) == (4,1)
@test num_allocs(cartpole) == 0

# Double Integrator
dim = 3
di = RobotZoo.DoubleIntegrator(dim)
n,m = size(di)
@test (n,m) == (6,3)
x,u = rand(di)
@test num_allocs(di) == 0

# Pendulum
pend = RobotZoo.Pendulum()
@test size(pend) == (2,1)
@test num_allocs(pend) == 0

# Planar Rocket

# Quadrotor
quad = RobotZoo.Quadrotor()
@test size(quad) == (13,4)
@test num_allocs(quad) == 0

# Yak Plane
yak = RobotZoo.YakPlane(MRP{Float64})
@test size(yak) == (12,4)
@test num_allocs(yak) == 0

# Test other functions
dt = 0.1
car = RobotZoo.DubinsCar()
n,m = size(car)
@test zeros(car) == (zeros(n), zeros(m))
@test zeros(Int,car)[1] isa SVector{n,Int}
@test fill(car,0.1) == (fill(0.1,n), fill(0.1,m))
@test ones(Float32,car)[2] isa SVector{m,Float32}

# Test default integrator
x,u = rand(car)
z = KnotPoint(x,u,dt)
@test discrete_dynamics(car, z) == discrete_dynamics(RK3, car, z)


text:
@testset "core.jl" begin

@testset "Trainer" begin
    _trainer = Trainer(;
        policy=BasicDQNLearner(
            approximator = NeuralNetworkApproximator(
                model = Chain(
                    Dense(4, 128, relu; initW = glorot_uniform),
                    Dense(128, 128, relu; initW = glorot_uniform),
                    Dense(128, 2; initW = glorot_uniform),
                ) |> cpu,
                optimizer = ADAM(),
            ),
            loss_func = huber_loss,
        )
    )

    trainer = actor(_trainer)

    put!(trainer, FetchParamMsg())
    ps = take!(self())
    original_sum = sum(sum, ps.data)

    for x in ps.data
        fill!(x, 0.)
    end

    put!(trainer, FetchParamMsg())
    ps = take!(self())
    new_sum = sum(sum, ps.data)

    # make sure no state sharing between messages
    @test original_sum == new_sum

    batch_data = (
        state = rand(4, 32),
        action = rand(1:2, 32),
        reward = rand(32),
        terminal = rand(Bool, 32),
        next_state = rand(4,32),
        next_action = rand(1:2, 32)
    )

    put!(trainer, BatchDataMsg(batch_data))

    put!(trainer, FetchParamMsg())
    ps = take!(self())
    updated_sum = sum(sum, ps.data)
    @test original_sum != updated_sum
end

@testset "TrajectoryManager" begin
    _trajectory_proxy = TrajectoryManager(
        trajectory = CircularSARTSATrajectory(;capacity=5, state_type=Any, ),
        sampler = UniformBatchSampler(3),
        inserter = NStepInserter(),
    )

    trajectory_proxy = actor(_trajectory_proxy)

    # 1. init traj for testing
    traj = CircularCompactSARTSATrajectory(
        capacity = 2,
        state_type = Float32,
        state_size = (4,),
    )
    push!(traj;state=rand(Float32, 4), action=rand(1:2))
    push!(traj;reward=rand(), terminal=rand(Bool),state=rand(Float32, 4), action=rand(1:2))
    push!(traj;reward=rand(), terminal=rand(Bool),state=rand(Float32, 4), action=rand(1:2))

    # 2. insert
    put!(trajectory_proxy, InsertTrajectoryMsg(deepcopy(traj)))  #!!! we used deepcopy here

    # 3. make sure the above message is already been handled
    put!(trajectory_proxy, PingMsg())
    take!(self())

    # 4. test that updating traj will not affect data in trajectory_proxy
    s_tp = _trajectory_proxy.trajectory[:state]
    s_traj = traj[:state]

    @test s_tp[1] == s_traj[:, 1]

    push!(traj;reward=rand(), terminal=rand(Bool),state=rand(Float32, 4), action=rand(1:2))

    @test s_tp[1] != s_traj[:, 1]

    s = sample(_trajectory_proxy.trajectory, _trajectory_proxy.sampler)
    fill!(s[:state], 0.)
    @test any(x -> sum(x) == 0, s_tp) == false  # make sure sample create an independent copy
end

@testset "Worker" begin
    _worker = Worker() do worker_proxy
        Experiment(
            Agent(
                policy = StaticPolicy(
                        QBasedPolicy(
                        learner = BasicDQNLearner(
                            approximator = NeuralNetworkApproximator(
                                model = Chain(
                                    Dense(4, 128, relu; initW = glorot_uniform),
                                    Dense(128, 128, relu; initW = glorot_uniform),
                                    Dense(128, 2; initW = glorot_uniform),
                                ) |> cpu,
                                optimizer = ADAM(),
                            ),
                            loss_func = huber_loss,
                        ),
                        explorer = EpsilonGreedyExplorer(
                            kind = :exp,
                            ϵ_stable = 0.01,
                            decay_steps = 500,
                        ),
                    ),
                ),
                trajectory = CircularCompactSARTSATrajectory(
                    capacity = 10,
                    state_type = Float32,
                    state_size = (4,),
                ),
            ),
            CartPoleEnv(; T = Float32),
            ComposedStopCondition(
                StopAfterStep(1_000),
                StopSignal(),
            ),
            ComposedHook(
                UploadTrajectoryEveryNStep(mailbox=worker_proxy, n=10, sealer=x -> InsertTrajectoryMsg(deepcopy(x))),
                LoadParamsHook(),
                TotalRewardPerEpisode(),
            ),
            "experimenting..."
        )
    end

    worker = actor(_worker)
    tmp_mailbox = Channel(100)
    put!(worker, StartMsg(tmp_mailbox))
end

@testset "WorkerProxy" begin
    target = RemoteChannel(() -> Channel(10))
    workers = [RemoteChannel(()->Channel(10)) for _ in 1:10]
    _wp = WorkerProxy(workers)
    wp = actor(_wp)

    put!(wp, StartMsg(target))
    for w in workers
        # @test take!(w).args[1] === wp
        @test Distributed.channel_from_id(remoteref_id(take!(w).args[1])) === Distributed.channel_from_id(remoteref_id(wp))
    end

    msg = InsertTrajectoryMsg(1)
    put!(wp, msg)
    @test take!(target) === msg

    for w in workers
        put!(wp, FetchParamMsg(w))
    end
    # @test take!(target).from === wp
    @test Distributed.channel_from_id(remoteref_id(take!(target).from)) === Distributed.channel_from_id(remoteref_id(wp))

    # make sure target only received one FetchParamMsg
    msg = PingMsg()
    put!(target, msg)
    @test take!(target) === msg

    msg = LoadParamMsg([])
    put!(wp, msg)
    for w in workers
        @test take!(w) === msg
    end
end

@testset "Orchestrator" begin
    # TODO
    # Add an integration test
end

end


text:
@register Base.getindex(x,i::Integer)
@register Base.getindex(x,i)
@register Base.binomial(n,k)

@register Base.signbit(x)
ModelingToolkit.derivative(::typeof(signbit), args::NTuple{1,Any}, ::Val{1}) = 0
ModelingToolkit.derivative(::typeof(abs), args::NTuple{1,Any}, ::Val{1}) = IfElse.ifelse(signbit(args[1]),-one(args[1]),one(args[1]))
function ModelingToolkit.derivative(::typeof(min), args::NTuple{2,Any}, ::Val{1})
    x, y = args
    IfElse.ifelse(x < y, one(x), zero(x))
end
function ModelingToolkit.derivative(::typeof(min), args::NTuple{2,Any}, ::Val{2})
    x, y = args
    IfElse.ifelse(x < y, zero(y), one(y))
end
function ModelingToolkit.derivative(::typeof(max), args::NTuple{2,Any}, ::Val{1})
    x, y = args
    IfElse.ifelse(x > y, one(x), zero(x))
end
function ModelingToolkit.derivative(::typeof(max), args::NTuple{2,Any}, ::Val{2})
    x, y = args
    IfElse.ifelse(x > y, zero(y), one(y))
end

@register IfElse.ifelse(x,y,z::Any)
ModelingToolkit.derivative(::typeof(IfElse.ifelse), args::NTuple{3,Any}, ::Val{1}) = 0
ModelingToolkit.derivative(::typeof(IfElse.ifelse), args::NTuple{3,Any}, ::Val{2}) = IfElse.ifelse(args[1],1,0)
ModelingToolkit.derivative(::typeof(IfElse.ifelse), args::NTuple{3,Any}, ::Val{3}) = IfElse.ifelse(args[1],0,1)

ModelingToolkit.@register Base.rand(x)
ModelingToolkit.@register Base.randn(x)

ModelingToolkit.@register Distributions.pdf(dist,x)
ModelingToolkit.@register Distributions.logpdf(dist,x)
ModelingToolkit.@register Distributions.cdf(dist,x)
ModelingToolkit.@register Distributions.logcdf(dist,x)
ModelingToolkit.@register Distributions.quantile(dist,x)

ModelingToolkit.@register Distributions.Uniform(mu,sigma)
ModelingToolkit.@register Distributions.Normal(mu,sigma)

@register ∈(x::Num, y::AbstractArray)
@register ∪(x, y)
@register ∩(x, y)
@register ∨(x, y)
@register ∧(x, y)
@register ⊆(x, y)

function LinearAlgebra.det(A::AbstractMatrix{<:Num}; laplace=true)
    if laplace
        n = LinearAlgebra.checksquare(A)
        if n == 1
            return A[1, 1]
        elseif n == 2
            return A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1]
        else
            temp = 0
            # Laplace expansion along the first column
            M′ = A[:, 2:end]
            for i in axes(A, 1)
                M = M′[(1:n) .!= i, :]
                d′ = A[i, 1] * det(M)
                if iseven(i)
                    temp = iszero(temp) ? d′ : temp - d′
                else
                    temp = iszero(temp) ? d′ : temp + d′
                end
            end
        end
        return temp
    else
        if istriu(A) || istril(A)
            return det(UpperTriangular(A))
        end
        return det(lu(A; check = false))
    end
end


text:
module Giggle

#PACKAGES USED
using Combinatorics
using DataFrames
using Distances
using Gurobi
using JuMP
using MathOptInterface
using Random
using UUIDs
using XLSX
using Query

#STRUCTS USED
struct vtx
    #IDENTIFIERS
    name::String
    type::String

    #PARAMS
    x::Float64 #x coor
    y::Float64 #y coor
    MAX::Float64 #max inventory level
    MIN::Float64 #min inventory level
    START::Float64 #starting inventory level

    #COSTS
    h::Float64
end

struct veh
    #IDENTIFIERS
    name::String
    type::String

    #PARAMS
    cover::Vector{Int64}
    Q::Int64
    start::Int64
    freq::Int64

    #COSTS
    varq::Float64
    vardq::Float64
    vard::Float64
    fix::Float64
end

struct col
    #QUANTITY RELATED
    q::JuMP.Containers.DenseAxisArray
    u::JuMP.Containers.DenseAxisArray
    v::JuMP.Containers.DenseAxisArray

    #0-1 VARS
    p::JuMP.Containers.DenseAxisArray
    y::JuMP.Containers.DenseAxisArray
    z::JuMP.Containers.DenseAxisArray
end

struct dval
    λ::JuMP.Containers.DenseAxisArray
    δ::JuMP.Containers.DenseAxisArray
end

struct stabilizer
    slackCoeff::Float64
    slackLim::JuMP.Containers.DenseAxisArray
    surpCoeff::Float64
    surpLim::JuMP.Containers.DenseAxisArray
end

struct bound
    idx::NamedTuple
    val::Int64
end

struct node
    #IDENTIFIERS
    parent::UUID
    self::UUID

    #BASE
    base::NamedTuple
    bounds::Vector{bound}

    #DYNAMIC SETS
    columns::Vector{col}

    #SUPPORT STRUCTURE
    stblzr::stabilizer
    status::Vector{String}
end

#CONSTANTS USED THROUGHOUT
const M = 9999999
const rng = MersenneTwister(1234)

#function to generate base data
function base(path::String)
    #READ WORKSHEET
    xf = XLSX.readxlsx(path)

    #DATAFRAME DICT
    data = Dict{Symbol,DataFrame}()

    #TURN SHEETS INTO DATAFRAME
    for sheets in XLSX.sheetnames(xf)
        #TRANSFORM TABLE INTO DATAFRAME
        df = DataFrame(XLSX.gettable(xf[sheets])...)

        #DEFINE THE NAME FROM THE WORKSHEET
        data[Symbol(sheets)] = df
    end

    #PROCESS VERTICES
    V = Dict{Int64,vtx}()

    #iterate over  each row of vertices
    for v in eachrow(data[:vertices])
        #direct input intu struct vtx
        V[v.id] = vtx(
            v.name,v.type,
            v.x,v.y,
            v.MAX,v.MIN,v.START,
            v.h
        )
    end

    #PROCESS VEHICLES
    K = Dict{Int64,veh}()

    #INITIATE INDEXING FOR VEHICLES
    idx = 0

    #iterate over each row of vehicles
    for v in eachrow(data[:vehicles])
        #convert string into vector{Int64}
        v.cover = parse.(Int64,split(v.cover))
        v.loadp = parse.(Int64,split(v.loadp))

        for f in v.loadp
            idx += 1
            K[idx] = veh(
                v.name,v.type,
                v.cover,v.Q,f,v.freq,
                v.varq,v.vardq,v.vard,v.fix
            )
        end
    end

    #EXTRACT T SET
    t_start = data[:periods].start[end] #STARTING MONTH
    t_end = data[:periods].start[end] + data[:periods].T[end] - 1 #FINAL MONTH (CALCULATE BASED ON DURATION)
    T = [i for i in t_start:t_end]

    #EXTRACT DEMANDS
    d = JuMP.Containers.DenseAxisArray{Float64}(undef,collect(keys(V)),T)
    for i in keys(V)
        row = @from x in data[:demands] begin
            @where x.point == i
            @select x
            @collect DataFrame
        end

        for t in T
            d[i,t] = row[2:end-3][t][end]
        end
    end

    #GENERATE DISTANCE MATRIX (FROM V)
    earth = 6378.137
    dist = JuMP.Containers.DenseAxisArray{Float64}(undef,collect(keys(V)),collect(keys(V)))
    for i in keys(V), j in keys(V)
        if i != j
            dist[i,j] = haversine([V[i].x,V[i].y],[V[j].x,V[j].y],earth)
        else
            dist[i,j] = M
        end
    end

    #GENERATE G MATRIX
    G = JuMP.Containers.DenseAxisArray{Float64}(undef,collect(keys(V)),collect(keys(K)))
    G .= M
    for k in keys(K)
        seed = K[k].cover[findmin([dist[K[k].start,j] for j in K[k].cover])[2]]
        for x in K[k].cover
            if x != seed
                G[x,k] = (min(dist[K[k].start,x]+dist[x,seed]+dist[seed,K[k].start] , dist[K[k].start,seed]+dist[seed,x]+dist[x,K[k].start]) - (dist[K[k].start,seed] + dist[seed,K[k].start]))
            else
                G[x,k] = 0
            end
        end
    end

    #GENERATE DELIVERY COST MATRIX
    deli = JuMP.Containers.DenseAxisArray{Float64}(undef,collect(keys(V)),collect(keys(K)))
    deli .= M
    for k in keys(K)
        for i in K[k].cover
            deli[i,k] = K[k].varq + K[k].vardq * dist[K[k].start,i]
        end
    end

    #DATA GENERATION STATUS
    trucks = sort(collect(keys(filter(p -> last(p).type == "truck",K))))
    ships = sort(collect(keys(filter(p -> last(p).type == "ship",K))))
    trains = sort(collect(keys(filter(p -> last(p).type == "train",K))))
    println("there are $(length(K)) vehicle data points with the composition:")
    println("$(length(trucks)) truck(s) for index $(first(trucks)) to $(last(trucks))")
    println("$(length(ships)) ship(s) for index $(first(ships)) to $(last(ships))")
    println("$(length(trains)) train(s) for index $(first(trains)) to $(last(trains))")
    println()
    println("there are $(length(V)) vertex data points with the composition:")
    println("$(length(collect(keys(filter(p -> last(p).type == "source",V))))) source(s)")
    println("$(length(collect(keys(filter(p -> last(p).type == "point",V))))) point(s)")

    return (K=K,V=V,T=T,d=d,dist=dist,G=G,deli=deli)
end

#function to initialize stabilizer
function initStab(dt::NamedTuple,vSlCoeff::Float64,vSurpCoeff::Float64)
    slackCoeff = vSlCoeff
    surpCoeff = vSurpCoeff
    slackLim = JuMP.Containers.DenseAxisArray{Float64}(undef,keys(dt.V),dt.T)
    surpLim = JuMP.Containers.DenseAxisArray{Float64}(undef,keys(dt.V),dt.T)
    for i in keys(dt.V),t in dt.T
        slackLim[i,t] = abs(dt.d[i,t])
        surpLim[i,t] = abs(dt.d[i,t])
    end

    return stabilizer(slackCoeff,slackLim,surpCoeff,surpLim)
end

#function to create root node
function root(dt::NamedTuple;slCoeff::Float64,suCoeff::Float64)
    id = uuid1(rng)
    root = node(
        id,id,
        dt,Vector{bound}(),
        Vector{col}(),
        initStab(dt,slCoeff,suCoeff),
        ["UNVISITED"]
    )

    return root
end

#function to update stabilizer
function updateStab(stabilizers::stabilizer,param::Float64)
    for i in first(stabilizers.slackLim.axes),t in last(stabilizers.slackLim.axes)
        stabilizers.slackLim[i,t] = param * stabilizers.slackLim[i,t]
        if stabilizers.slackLim[i,t] < 1
            stabilizers.slackLim[i,t] = 0
        end
    end

    for i in first(stabilizers.surpLim.axes),t in last(stabilizers.surpLim.axes)
        stabilizers.surpLim[i,t] = param * stabilizers.surpLim[i,t]
        if stabilizers.surpLim[i,t] < 1
            stabilizers.surpLim[i,t] = 0
        end
    end

    return stabilizers
end

#===============MASTER FUNCTION SET=====#
function master(n::node;silent::Bool,env::Gurobi.Env)
    #build and bound
    mp = buildMaster(n;silent=silent,env=env)

    #solve the mp
    optimize!(mp)

    return mp
end

function buildMaster(n::node;silent::Bool,env::Gurobi.Env)
    #COLUMN LABELING
    R = Dict(1:length(n.columns) .=> n.columns)

    #MODEL DECLARATION
    mp = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(env)))
    if silent
        set_silent(mp)
    end

    #VARIABLE DECLARATION
    θ = @variable(mp, θ[keys(R),keys(n.base.K),n.base.T] >= 0) #LABEL 20
    I = @variable(mp, I[keys(n.base.V),vcat(first(n.base.T)-1,n.base.T)]) #NO BOUNDS (bounded below)

    #SLACK SURPLUS DECLARATION
    @variable(mp, 0 <= slack_I[i=keys(n.base.V),t=n.base.T] <= n.stblzr.slackLim[i,t])
    @variable(mp, 0 <= surp_I[i=keys(n.base.V),t=n.base.T] <= n.stblzr.surpLim[i,t])

    #CONSTRAINT 17 & 18 + STARTING INVENTORY
    @constraint(
        mp, λ[i=keys(n.base.V),t=n.base.T],
        I[i,t-1] + sum(R[r].q[i,k,t] * θ[r,k,t] for r in keys(R),k in keys(n.base.K)) + slack_I[i,t] - surp_I[i,t] == n.base.d[i,t] + I[i,t]
    )
    @constraint(mp, [i=keys(n.base.V),t=n.base.T], n.base.V[i].MIN <= I[i,t] <= n.base.V[i].MAX)
    @constraint(mp, [i=keys(n.base.V)], I[i,first(n.base.T)-1] == n.base.V[i].START)

    #CONVEXITY CONSTRAINT LABEL 19
    @constraint(mp, δ[k=keys(n.base.K),t=n.base.T], sum(θ[r,k,t] for r in keys(R)) <= n.base.K[k].freq)

    #BOUND GENERATOR
    #empty

    #OBJECTIVE FUNCTION
    begin
        @objective(mp, Min,
            sum(
                (
                    n.base.K[k].vard * g(R[r].y[:,k,t];n=n,k=k) +
                    sum(n.base.deli[i,k] * R[r].u[i,k,t] for i in n.base.K[k].cover) +
                    n.base.K[k].fix * R[r].z[n.base.K[k].start,k,t]
                ) * θ[r,k,t] for r in keys(R),k in keys(n.base.K),t in n.base.T
            ) +
            sum(n.base.V[i].h * I[i,t] for i in keys(n.base.V),t in n.base.T) +
            sum(n.stblzr.slackCoeff * slack_I[i,t] for i in keys(n.base.V),t in n.base.T) -
            sum(n.stblzr.surpCoeff * surp_I[i,t] for i in keys(n.base.V),t in n.base.T)
        )
    end

    return mp
end

function getDual(mp::Model)
    return duals = dval(dual.(mp.obj_dict[:λ]),dual.(mp.obj_dict[:δ]))
end
#===============MASTER FUNCTION SET=====#

#===============SUB FUNCTION SET=====#
function sub(n::node,duals::dval;silent::Bool,env::Gurobi.Env)
    #build and bound
    sp = buildSub(n,duals;silent=silent,env=env)

    #solve the sp
    optimize!(sp)

    return sp
end

function buildSub(n::node,duals::dval;silent::Bool,env::Gurobi.Env)
    #MODEL DECLARATION
    sp = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(env)))
    if silent
        set_silent(sp)
    end

    #VARIABLE DECLARATION
    q = @variable(sp, q[keys(n.base.V),keys(n.base.K),n.base.T],Int)
    u = @variable(sp, u[keys(n.base.V),keys(n.base.K),n.base.T] >= 0,Int)
    v = @variable(sp, v[keys(n.base.V),keys(n.base.K),n.base.T] >= 0,Int)
    y = @variable(sp, y[keys(n.base.V),keys(n.base.K),n.base.T], Bin)
    z = @variable(sp, z[keys(n.base.V),keys(n.base.K),n.base.T], Bin)
    p = @variable(sp, p[keys(n.base.V),keys(n.base.K),n.base.T], Bin)

    #CONSTRAINTS
    @constraint(sp, [k=keys(n.base.K),t=n.base.T], sum(q[i,k,t] for i in n.base.K[k].cover) == 0)
    @constraint(sp, [k=keys(n.base.K),i=n.base.K[k].cover,t=n.base.T], q[i,k,t] == u[i,k,t] - v[i,k,t])
    @constraint(sp, [k=keys(n.base.K),i=n.base.K[k].cover,t=n.base.T], u[i,k,t] <= n.base.K[k].Q * y[i,k,t])
    @constraint(sp, [k=keys(n.base.K),i=n.base.K[k].cover,t=n.base.T], v[i,k,t] <= n.base.K[k].Q * z[i,k,t])
    @constraint( #everything outside the possible starting point is zero
        sp, [k=keys(n.base.K),t=n.base.T],
        sum(z[i,k,t] for i in filter(x -> x != n.base.K[k].start,n.base.K[k].cover)) == 0
    )
    @constraint( #the starting point can be zero or one
        sp, [k=keys(n.base.K),t=n.base.T],
        z[n.base.K[k].start,k,t] <= 1
    )
    @constraint(sp, [k=keys(n.base.K),i=n.base.K[k].cover,t=n.base.T], p[i,k,t] == y[i,k,t] + z[i,k,t])

    #BOUND GENERATOR SUBPROBLEM
    #empty

    #OBJECTIVE FUNCTION
    begin
        @objective(sp,Min,
            sum(
                sum(n.base.K[k].vard * n.base.G[i,k] * y[i,k,t] for i in n.base.K[k].cover) +
                sum(n.base.deli[i,k] * u[i,k,t] for i in n.base.K[k].cover) +
                n.base.K[k].fix * z[n.base.K[k].start,k,t] for k in keys(n.base.K),t in n.base.T
            ) -
            sum(sum(q[i,k,t] * duals.λ[i,t] for i in n.base.K[k].cover) for k in keys(n.base.K),t in n.base.T) -
            sum(duals.δ[k,t] for k in keys(n.base.K),t in n.base.T)
        )
    end

    return sp
end

function getCol(sp::Model)
    #EXTRACT VARIABLES
    q = value.(sp.obj_dict[:q])
    u = value.(sp.obj_dict[:u])
    v = value.(sp.obj_dict[:v])
    p = value.(sp.obj_dict[:p])
    y = value.(sp.obj_dict[:y])
    z = value.(sp.obj_dict[:z])

    return col(q,u,v,p,y,z)
end

function realPrice(n::node,duals::dval;column::col)
    #CALCULATE REAL PRICING
    begin
        price = (
            sum(
                n.base.K[k].vard * g(column.y[:,k,t];n=n,k=k) +
                sum(n.base.deli[i,k] * column.u[i,k,t] for i in n.base.K[k].cover) +
                n.base.K[k].fix * column.z[n.base.K[k].start,k,t] for k in keys(n.base.K),t in n.base.T
            ) -
            sum(sum(column.q[i,k,t] * duals.λ[i,t] for i in n.base.K[k].cover) for k in keys(n.base.K),t in n.base.T) -
            sum(duals.δ[k,t] for k in keys(n.base.K),t in n.base.T)
        )
    end

    return price
end
#===============SUB FUNCTION SET=====#

#===============SUPPORT ESTIMATOR FOR TSP===#
function g(y::JuMP.Containers.DenseAxisArray;n::node,k::Int64)
    x = twoOpt([i for i in first(y.axes) if y[i] == 1.0];n=n,k=k)
    return x
end #CLEARED

function twoOpt(p::Vector{Int64};n::node,k::Int64) #INPUT IS A COLLECTION OF POINTS [15 22 23 24] ORDER IS ARBITRARY
    if length(p) > 0
        #DETERMINE START AND NODE
        start = n.base.K[k].start #passing k stop di sini
        nodes = vcat(start,p)

        #DO NEAREST NEIGHBORHOOD construction
        nnPath = nn(start,nodes;n=n)

        #RETURN ITS IMPROVEMENT
        return trav(improve(nnPath;n=n);n=n)
    else
        return 0
    end
end #CLEARED

function nn(start::Int64,nodes::Vector{Int64};n::node)
    tour = [start]

    unvisited = deepcopy(nodes)
    current_index = 1
    terminate = false

    while !terminate
        current_position = tour[current_index]

        #REMOVE CURRENT NODES AFTER FIRST MOVE
        to_remove = findfirst(x -> x==current_position,unvisited)
        splice!(unvisited,to_remove)

        #if no more unvisited, terminate
        if length(unvisited) == 0
            push!(tour,first(tour))
            terminate = true
        else
            terpilih = selectNn(current_position,unvisited;n=n)
            push!(tour,terpilih)
        end

        #UPDATE INDEX
        current_index += 1
    end

    return tour
end #CLEARED

function selectNn(current::Int64,unvisited::Vector{Int64};n::node)
    chosen = unvisited[1]

    for j in unvisited
        if n.base.dist[current,chosen] > n.base.dist[current,j]
            chosen = j #tuker kalo jaraknya lebih kecil
        end
    end

    return chosen
end #CLEARED

function improve(p::Vector{Int64};n::node) #INPUT: NNPATH i.e. [15 27 28 29 15]
    a = view(p,2:length(p)-1,1)

    for i in 1:length(a) - 1
        for k in i+1:length(a) - 1
            best = trav(p;n=n)
            reverse!(view(a,i:k,1))
            if trav(p;n=n) >= best
                reverse!(view(a,i:k,1))
            end
        end
    end

    return p
end #CLEARED

function trav(p::Vector{Int64};n::node) #INPUT: NNPATH i.e. [15 27 28 29 15]
    return sum(n.base.dist[p[i],p[i+1]] for i in 1:length(p)-1)
end
#===============SUPPORT ESTIMATOR FOR TSP===#

function check(model::Model)
    return value(sum(model.obj_dict[:slack_I])) + value(sum(model.obj_dict[:surp_I]))
end #CLEARED

function colGen(n::node,maxCG::Float64;silent::Bool,env::Gurobi.Env,track::Bool)
    #INITIALIZE
    terminate = false
    iter = 0

    while !terminate
        #GENERATE MASTER PROBLEM
        mp = master(n;silent=silent,env=env)

        #NODE PROCESSING STATUS CHECK
        if has_values(mp) && has_duals(mp)
            if track #tracking status
                println("obj: $(objective_value(mp))")
            end

            #EXTRACT DUAL
            duals = getDual(mp)

            #GENERATE SUBPROBLEM
            sp = sub(n,duals;silent=silent,env=env)

            #EXTRACT COLS & PRICE
            cols = getCol(sp)
            price = realPrice(n,duals;column=cols)
            slacksurp = check(mp)

            if track #tracking status
                println(price)
                println("nilai slack surp $slacksurp")
            end

            #NEGATIVE COLUMN CHECK
            if isapprox(price, 0, atol = 1e-8) || price > 0
                if isapprox(slacksurp, 0, atol = 1e-8)
                    terminate = true
                    push!(n.status,"EVALUATED")
                    if track
                        println("EVALUATED")
                    end
                else
                    updateStab(n.stblzr,0.4)
                    push!(n.status,"STABILIZED")
                    if track
                        println("STABILIZED")
                    end
                end
            else
                push!(n.columns,cols)
                push!(n.status,"ADD_COLUMN")
                if track
                    println("ADD_COLUMN")
                end
            end

            #UPDATE ITER
            iter += 1

            #MAXIMUM ITER
            if iter >= maxCG
                terminate = true
                pop!(n.columns)
                push!(n.status,"EVALUATED-TIME OUT")
                if track
                    println("EVALUATED-TIME OUT")
                end
            end
        else
            terminate = true
            push!(n.status,"NO_SOLUTION")
            if track
                println("NO_SOLUTION")
            end
        end
    end

    if n.status[end] == "EVALUATED" || n.status[end] == "EVALUATED-TIME OUT"
        println("NODE $(n.self) FINISHED.")
    else
        println("NODE $(n.self) FAILED.")
    end
end

export base
export root

export master
export buildMaster
export getDual

export sub
export buildSub
export getCol

export realPrice

export colGen

end


text:
using DataFrames
using Symbolics
using Interpolations
using DiffRules
using SymbolicUtils

# differentiate a function specified taken from tabular data
# data is an interpolated AbstractArray and the original data points
function dataDerivative( data::Tuple{AbstractArray,Vector{Float64}} )::Tuple{AbstractArray,Vector{Float64}}
    xs = data[ 2 ]
    ys = map( data[ 1 ], xs )

    if( length(xs) != length(ys))
        throw("Data doesn't have same number of xs and ys.")
    end

    l = length( ys )
    nys = if l > 2
        vcat(
            [ (ys[2] - ys[1])/(xs[2] - xs[1]) ],
            [ (ys[i+1] - ys[i-1])/(xs[i+1] - xs[i-1]) for i in 2:(l-1) ],
            [ (ys[l] - ys[l-1])/(xs[l] - xs[l-1])]
        )
    elseif l == 2
        [
            (ys[2] - ys[1])/(xs[2] - xs[1]),
            (ys[2] - ys[1])/(xs[2] - xs[1])
        ]
    else
        throw( "Data must have atleast two points to be able to take derivative." )
    end

    nData = ( LinearInterpolation( xs, nys; extrapolation_bc = Line() ), xs )
    return( nData )
end

# just return the value at index x from the abstractvector in data
function dataFunction( data::Tuple{AbstractArray,Vector{Float64}}, x::Float64 )::Float64
    data[ 1 ](x)
end
# if this derivative is not specified correctly then there may be a weird error when trying to expand_derivatives on this function
# Symbolics.derivative( ::typeof(dataFunction), args::NTuple{2,Any}, ::Val{2} ) = SymbolicUtils.Term( dataFunction, [ dataDerivative(args[1]), args[2] ] )
# Symbolics.@register dataFunction( data::Tuple{AbstractArray,Vector{Float64}}, x )


text:
# clones a table
# If the flush option is not enabled, then any data the source table currently has in memory will not exist in the clone.
# A cloned table copies the configuration of the source table.
# However the permissions of the source table are not copied to the clone.
# After a clone is created, only the user that created the clone can read and write to it.
function table_clone(session::AccumuloSession, tablename::String, newtablename::String; flush::Bool=true, set_properties::Dict=Dict(), exclude_properties::SET=())
    setprops = Dict{String,String}(string(n)=>string(v) for (n,v) in set_properties)
    exclprops = Set(String[string(x) for x in exclude_properties])
    cloneTable(client(session), handle(session), tablename, newtablename, flush, setprops, exclprops)
end

# create, delete, rename tables
#
# When a table is created, by default it is configured to use the VersioningIterator and keep one version. A table can be created without the VersioningIterator.
# Accumulo can be configured to return the top k versions, or versions later than a given date. The default is to return the one most recent version.
# The version policy can be changed by changing the VersioningIterator options with the `config` api.
#
# Using logical timestamps ensures that timestamps set by Accumulo always move forward. This helps avoid problems caused by TabletServers that have different time settings.
# The default is to use milliseconds as timestamp.
function table_create(session::AccumuloSession, tablename::String; versioning::Bool=true, logical_timestamp::Bool=false)
    time_type = logical_timestamp ? TimeType.LOGICAL : TimeType.MILLIS
    createTable(client(session), handle(session), tablename, versioning, time_type)
end
table_delete(session::AccumuloSession, tablename::String) = deleteTable(client(session), handle(session), tablename)
table_rename(session::AccumuloSession, tablename::String, new_name::String) = renameTable(client(session), handle(session), tablename, new_name)

# prints how much space, in bytes, is used by files referenced by a table.  When multiple tables are specified it prints how much space, in bytes, is used by files shared between tables, if any.
function table_du(session::AccumuloSession, tables::String...)
    tableset = Set(String[string(x) for x in tables])
    usage = getDiskUsage(client(session), handle(session), tableset)
    (usage[1]).usage
end

table_exists(session::AccumuloSession, table::String) = tableExists(client(session), handle(session), table)

# get or set table specific properties
table_config(session::AccumuloSession, tablename::String) = getTableProperties(client(session), handle(session), tablename)
table_config!(session::AccumuloSession, tablename::String, property::String, value::String) = setTableProperty(client(session), handle(session), tablename, property, value)

# limit versions in a table
function table_versions!(session::AccumuloSession, tablename::String, nversions::Integer)
    (nversions > 0) || error("at least 1 version required")
    viter = iter("vers", "org.apache.accumulo.core.iterators.user.VersioningIterator", 10, Dict("maxVersions"=>nversions))
    add_iter(session, tablename, viter, (IteratorScope.MINC, IteratorScope.MAJC, IteratorScope.SCAN))
end

# import/export tables
#
# Accumulo supports exporting tables for the purpose of copying tables to another cluster.
# Exporting and importing tables preserves the tables configuration, splits, and logical time.
# Tables are exported and then copied via the hadoop distcp command.
# To export a table, it must be offline and stay offline while discp runs. The reason it needs to stay offline is to prevent files from being deleted.
# A table can be cloned and the clone taken offline inorder to avoid losing access to the table.
table_export(session::AccumuloSession, tablename::String, export_dir::String) = exportTable(client(session), handle(session), tablename, export_dir)
table_import(session::AccumuloSession, tablename::String, import_dir::String) = importTable(client(session), handle(session), tablename, import_dir)

# take tables online/offline
table_offline(session::AccumuloSession, tablename::String; wait::Bool=true) = offlineTable(client(session), handle(session), tablename, wait)
table_online(session::AccumuloSession, tablename::String; wait::Bool=true) = onlineTable(client(session), handle(session), tablename, wait)

# get a list of all existing tables
tables(session::AccumuloSession) = listTables(client(session), handle(session))

#=
namespaces not supported yet
# displays a list of all existing namespaces
function namespaces()
end

# creates a new namespace
function createnamespace()
end

# deletes a namespace
function deletenamespace()
end

# renames a namespace
function renamenamespace()
end
=#


text:
using Test
using LinearAlgebra

CASE = joinpath(dirname(@__FILE__), "..", "data", "case9")

@testset "OPFData" begin
    data = ExaTron.opf_loaddata(CASE)
    @test isa(data, ExaTron.OPFData)
end

@testset "AdmmEnv" begin
    use_gpu = false
    rho_pq = 1.0
    rho_va = 1.0
    env = ExaTron.AdmmEnv(CASE, use_gpu, rho_pq, rho_va)

    @test env.case == CASE
    @test env.solution.status == ExaTron.INITIAL

    @testset "AdmmEnv.Model" begin
        model = env.model
        gen_model = model.gen_mod
        ngen = gen_model.ngen
        @test length(gen_model.pgmin) == length(gen_model.pgmax) == length(gen_model.qgmin) == length(gen_model.qgmax) == ngen
        @test length(gen_model.c2) == length(gen_model.c1) == length(gen_model.c0) == ngen

        nlines = model.nline
        @test length(model.FrBound) == length(model.ToBound) == 2 * nlines
        @test length(model.YshR) == length(model.YshI) == length(model.YffR) ==
              length(model.YffI) == length(model.YftR) == length(model.YftI) ==
              length(model.YttR) == length(model.YttI) == length(model.YtfR) ==
              length(model.YtfI) == nlines

        nbus = model.nbus
        @test length(model.Pd) == length(model.Qd) == nbus
    end
end

@testset "One-level ADMM algorithm" begin
    # NB: Need to run almost 2,000 iterations to reach convergence with this
    # set of parameters.
    env = ExaTron.admm_rect_gpu(CASE; verbose=0, iterlim=2000, rho_pq=400.0, rho_va=40000.0)
    @test isa(env, ExaTron.AdmmEnv)

    model = env.model
    ngen = model.gen_mod.ngen

    par = env.params
    sol = env.solution

    # Check results
    pg = ExaTron.active_power_generation(env)
    qg = ExaTron.reactive_power_generation(env)

    @test sol.status == ExaTron.HAS_CONVERGED
    # Test with solution returned by PowerModels + Ipopt
    @test sol.objval ≈ 5296.6862 rtol=1e-4
    @test pg ≈ [0.897987, 1.34321, 0.941874] rtol=1e-4
    @test qg ≈ [0.1296564, 0.00031842, -0.226342] rtol=1e-2

    # Test restart API
    ExaTron.admm_restart!(env)
end

@testset "Two-level ADMM algorithm" begin
    # Two-level algorithm has around two times more coupling constraints
    # than the one-level algorithm, so it is expected to be less accurate
    # than the one level algorithm.
    env = ExaTron.admm_rect_gpu_two_level(
        CASE;
        verbose=0, rho_pq=1000.0, rho_va=1000.0, inner_iterlim=2000, outer_iterlim=1500, outer_eps=1e-6
    )
    @test isa(env, ExaTron.AdmmEnv)

    model = env.model
    ngen = model.gen_mod.ngen

    par = env.params
    sol = env.solution

    # Check results
    pg = ExaTron.active_power_generation(env)
    qg = ExaTron.reactive_power_generation(env)

    @test sol.status == ExaTron.HAS_CONVERGED
    # # Test with solution returned by PowerModels + Ipopt
    @test_broken sol.objval ≈ 5296.6862 rtol=1e-4
    @test_broken pg ≈ [0.897987, 1.34321, 0.941874] rtol=1e-4
    @test qg ≈ [0.1296564, 0.00031842, -0.226342] atol=1e-2

    # Test restart API
    ExaTron.admm_restart!(env)
end


text:
module AtomicPotentials
using AtomicLevels

abstract type AbstractPotential{T} end

# * Table of elements

include("table_of_elements.jl")

function element_number(s::Symbol)
    i = findfirst(e -> e[1] == s, table_of_elements)
    i === nothing &&
        throw(ArgumentError("Invalid element $s"))
    i
end

# * Point charge nucleus

struct PointCharge{T} <: AbstractPotential{T}
    Z::T
end

PointCharge(s::Symbol) = PointCharge(element_number(s))

macro pc_str(s)
    :(PointCharge(Symbol($s)))
end

function Base.show(io::IO, p::PointCharge{I}) where {I<:Integer}
    element = table_of_elements[p.Z]
    write(io, "Z = $(p.Z) [$(element[2][1])]")
end

Base.show(io::IO, p::PointCharge) =
    write(io, "Z = $(p.Z)")

(p::PointCharge{T})(::O, r::U) where {T,O,U} = -p.Z/r
(p::PointCharge{T})(::O, r::VU) where {T,O,U,VU<:AbstractVector{U}} = -p.Z./r

charge(p::PointCharge) = p.Z

ground_state(p::PointCharge{<:Integer}) =
    table_of_elements[p.Z][2][2]

ground_state(p::PointCharge) =
    throw(ArgumentError("Ground state configuration for nuclear charge of Z = $(p.Z) unknown"))

islocal(::PointCharge) = true

# * Yukawa
struct Yukawa{T} <: AbstractPotential{T}
    g::T
    α::T
end

charge(::Yukawa) = 1
ground_state(::Yukawa) = c"1s"

(p::Yukawa{T})(::O, r::U) where {T,O,U} = -p.g^2*exp(-p.α*r)/r

islocal(::Yukawa) = true

# * Exports

export AbstractPotential, element_number, PointCharge, @pc_str, charge, ground_state, islocal, Yukawa

end # module


text:
promote_u0(u0::AbstractArray{<:ForwardDiff.Dual},p::AbstractArray{<:ForwardDiff.Dual},t0) = u0
promote_u0(u0,p::AbstractArray{<:ForwardDiff.Dual},t0) = eltype(p).(u0)
promote_u0(u0,p::NTuple{N,<:ForwardDiff.Dual},t0) where N = eltype(p).(u0)
promote_u0(u0,p::ForwardDiff.Dual,t0) where N = eltype(p).(u0)

function promote_tspan(u0::AbstractArray{<:ForwardDiff.Dual},p,tspan::Tuple{<:ForwardDiff.Dual,<:ForwardDiff.Dual},prob,kwargs)
  return tspan
end

function promote_tspan(u0::AbstractArray{<:ForwardDiff.Dual},p,tspan,prob,kwargs)
  if (haskey(kwargs,:callback) && has_continuous_callback(kwargs[:callback])) ||
     (haskey(prob.kwargs,:callback) && has_continuous_callback(prob.kwargs[:callback]))

    return eltype(u0).(tspan)
  else
    return tspan
  end
end

value(x::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N} = V
value(x::ForwardDiff.Dual) = value(ForwardDiff.value(x))

@inline fastpow(x::ForwardDiff.Dual, y::ForwardDiff.Dual) = x^y

sse(x::Number) = x^2
sse(x::ForwardDiff.Dual) = sse(ForwardDiff.value(x)) + sum(sse, ForwardDiff.partials(x))
totallength(x::Number) = 1
totallength(x::ForwardDiff.Dual) = totallength(ForwardDiff.value(x)) + sum(totallength, ForwardDiff.partials(x))
totallength(x::AbstractArray) = sum(totallength,x)

@inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual,::Any) = sqrt(sse(u))
@inline ODE_DEFAULT_NORM(u::AbstractArray{<:ForwardDiff.Dual},t::Any) = sqrt(sum(sse,u) / totallength(u))
@inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual,::ForwardDiff.Dual) = sqrt(sse(u))
@inline ODE_DEFAULT_NORM(u::AbstractArray{<:ForwardDiff.Dual},::ForwardDiff.Dual) = sqrt(sum(x->sse(x),u) / totallength(u))

if !hasmethod(nextfloat, Tuple{ForwardDiff.Dual})
  # Type piracy. Should upstream
  Base.nextfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(nextfloat(d.value), d.partials)
  Base.prevfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(prevfloat(d.value), d.partials)
end

# bisection(f, tup::Tuple{T,T}, t_forward::Bool) where {T<:ForwardDiff.Dual} = find_zero(f, tup, Roots.AlefeldPotraShi())


text:
using Mimi

# Load RICE+AIR source code.
include("src/mRICE2010.jl")

m = get_model()
run(m)

explore(m)

# say you want to change some parameter
#like the discount rate
m[:welfare,:rho] = 1.5
#or the mitigation fraction, this is how we input the control variable in the optimisation
m[:emissions,:MIU] = 0.5*ones(60,12)

max_stars_repo_name: fdennig/RICEupdate

max_issues_repo_name: fdennig/RICEupdate

max_forks_repo_name: fdennig/RICEupdate


text:
using ITensors
using ITensorVisualization
using Graphs

tn = itensornetwork(grid((3, 3, 3)))
edge_labels = (; dims=false)
@visualize fig tn ndims = 3 edge_labels = edge_labels vertex_size = 400

fig


text:
module Flow

module Args
	import ...Utilities
	import ...Utilities.Args: Arg
	#export accepted_arg_list, throw_if_incompatible_options
	# No extra flags for this model yet.
	Utilities.Args.accepted_arg_list(::Val{:Flow}) = Vector{Arg}()
	function Utilities.Args.throw_if_incompatible_options(
		::Val{:Flow}, p_args
	)
		nothing
	end
end

# Used for building the model, used by module Format below.
include("Enumeration.jl")
using .Enumeration

# Just include it to make it available to users of this module.
include("Format.jl")

using JuMP

import ..build_model
import ..BUILT_MODEL
"""
    build_model(::Val{:Flow}, model, d, p, l, w, L, W, [options])

The Flow model is a variant of the PPG2KP model using flow constraints.
The result is a considerably larger model with a slight enhancement
on the relaxation bound.

# Arguments

The required arguments are the same as all `build_model` methods, and
the optional argument `options` is always ignored.

# High-Level Model Description

Both length and width of the original plate may be seen as line segments.
From these line segments, let us just consider a subset of the points:
the first point (zero), the last point (L or W), and the ones in the middle
that match a normal cut position.

For each width cut position, there is an instance of the length point set
representing plates with all possible lengths but with a fixed width.
The analogue (each length cut position ... width point set) also exists.
Each of these point set instances is the vertex set of a strongly connected
directed graph.
There are no edges connecting one of such graphs with another (they interact by
special constrainsts that cannot be seen as 'flow constraints'/edges).

The variables of the model are the edges between those graphs.
The edges are may be divided into forward edges and backward/circulation
edges.
The forward edges go from a node representing a point closer to the start
of the line segment to a node representing a point farther from the start.
The backward/circulation edges always go from a non-first (non-zero) point to
the first/zero point of the graph.
The forward edges represent cutting part of the plate for use or trimming.
The backard edges represent the availibity of a plate of some dimensions
to be cut in some orientation.
For example, if there are two units of flow flowing through a forward edge that
goes from point length 15 to point length 35 of the set with fixed width 50,
then there are at least two plates of width 50 that have a subplate of length
20 (and width 50) cut from two parallel cuts, one in position 15 and another
in 35. The two parent plates may have the same dimensions (i.e., their length
is also equal to each other), or they may have different total length (they
only share the same width); for an example, one of them may stop at length 35
and the other may have length 1000. Considering this same example, if there are
a plate of size 35x50 and another of size 1000x50 available, then there is one
of flow in the backward edge from point length 35 to point zero and from point
length 1000 to point zero (all inside the set with fixed width 50).

Finally, forward edges may be subdivided into waste edges, piece edges and
subplate edges.
Waste edges do nothing besides letting the flow pass through them.
The variables representing piece edges are in the objective function with
coefficient equal to the profit of the corresponding piece.
Subplate edges "introduce" new flow in the system, they do not only allow the
flow to pass through them but also increase the flow in a correspondent
circulation edge of the (other) graph that represents the plate they are
extracting, but to be cut in the opposite orientation.
For example, an increase in the flow of the subplate edge from length 50 to
length 150 of the set with fixed width 200 will increase by the same amount the
flow of the backward edge from width 200 to width zero of the set with fixed
length 100.

The number and dimension of the original plates available is given by all
backward edges that have their value fixed to a non-zero value a priori.

# Slightly Misleading Low-Level Explanation

## Variables

* `edge[i]`: Positive Integer. The only variable. Measures the amount of flow
  passing through the edge. Has all edges indexed by their globalized index.
  Uses the many reverse indexes created by the enumeration process to build the
  model.

## Objective function:

Maximize the profit of the pieces cut.

`sum(p[CP] * edge[1:number_of_pieces])`

## Constraints

* There is exactly one of the original plate, so its circulation arc is
  available since start, and cannot ever go over one. It must be one.
  - `sum(edge[is the backward edge of the original plate]) == 1`
* The flow must go round.
  - `sum(edge[all edges ending in node x]) == sum(edge[all edges starting at
    node x]) : for every node x`
* If you cut a subplate anywhere then you allow its flow to circle again.
  - `sum(edge[a backward edge]) <= sum(edge[a forward edge that cuts the
    subplate corresponging to that backward edge]) : for each backward edges`
* The number of pieces of some type is always less than or equal to its demand.
  - `edge[i] <= d[i] : for each i in 1:number_of_pieces`

## Dummy Plate Identifiers:
* 0 - No plate. Waste and backward edges use this value in the `back` field.
* 1:n - A dummy edge corresponding to a piece index, used in the objective
  function. Edges cutting pieces do not actually appear in the objective
  function but, in fact, allow more flow to pass through by some of these
  dummy piece edges.
* n+1 - A dummy edge that is upper bounded at one and allows one extra flow
  to the original plate default cutting (vertical) backward edge (i.e.,
  n+2). Necessary to avoid n+2 to be restricted by a possible edge that
  cuts the whole original plate alternative cutting (horizontal). In other
  words, without this, if a piece has the same length as the plate, then
  the model has a deadlock problem, and will give a zero objective problem.
* n+2 - The original plate default cutting (vertical) backward edge.
* n+3 - The original plate vertical forward edge that takes the whole plate
  and allows it to be horizontally cut (i.e., allow flow to n+2).
* n+4 - The original plate alternative cutting (horizontal) backward edge.
  It is allowed by cutting n+2.
* n+5:m - The remaining intermediary plates.
"""
function build_model(
	::Val{:Flow}, model, d :: Vector{D}, p :: Vector{P},
	l :: Vector{S}, w :: Vector{S}, L :: S, W :: S,
	options :: Dict{String, Any} = Dict{String, Any}()
) where {D, S, P}
	# TODO: fix (L, W, l, w, d) = (20, 5, [20, 4, 18], [4, 5, 2], [1, 1, 1])
	@assert length(d) == length(l)
	@assert length(d) == length(w)
	@assert length(d) == length(p)
	num_piece_types = convert(D, length(d))

	N = E = P
	# TODO: check if N and E will be bubbled up.
	# NOTE: the profits are not needed, the first num_piece_types edges will
	# be dummy edges that have the profit value in the objective function.
	nodes_data, edges_data, last_gnode_idx, last_gedge_idx = gen_all_edges(
		N, E, d, l, w, L, W
	)
	#@show edges_data

	head2indxs = [Vector{E}() for _ = 1:last_gnode_idx]
	tail2indxs = [Vector{E}() for _ = 1:last_gnode_idx]
	back2indxs = [Vector{E}() for _ = 1:last_gedge_idx]
	for e in edges_data
		@assert e.indx <= last_gedge_idx
		@assert e.head <= last_gnode_idx
		iszero(e.head) && @show e
		@assert e.tail <= last_gnode_idx
		iszero(e.tail) && @show e
		push!(head2indxs[e.head], e.indx)
		push!(tail2indxs[e.tail], e.indx)
		!iszero(e.back) && push!(back2indxs[e.back], e.indx)
	end
	# Create the dummy edge that enables the "original plate default orientation
	# backward edge" to pass one unit of flow.
	dummy_vroot_enabler = num_piece_types + one(D)
	vroot_back_edge = dummy_vroot_enabler + one(D)
	#@assert isempty(back2indxs[vroot_back_edge])
	push!(back2indxs[vroot_back_edge], dummy_vroot_enabler)

	@variable(model, edge[1:last_gedge_idx] >= 0, Int)

	fitting = (pii for pii = 1:num_piece_types if (l[pii] <= L && w[pii] <= W))

	if !isempty(fitting)
		@objective(model, Max, sum(p[pii] * edge[pii] for pii in fitting))
	end

	# There is just one of the original plate. Sets the dummy edge to one and
	# consequently allow one unit of flow to the original plate default
	# orientation backward edge.
	@constraint(model, edge[dummy_vroot_enabler] == 1)

	for i = 1:last_gnode_idx
		# The flow model should not have vertexes that have no arc arriving
		# or leaving it (flow must be able to reach the vertex and leave it,
		# however in the rare case of an instance with pieces that do not fit
		# the plate, or a empty set of pieces, this may happen.
		#=
		if isempty(head2indxs[i]) || isempty(tail2indxs[i])
			@show i
			@show head2indxs[i]
			@show tail2indxs[i]
		end
		@assert !isempty(head2indxs[i]) && !isempty(tail2indxs[i])
		=#
		@constraint(model,
			sum(edge[head2indxs[i]]) == sum(edge[tail2indxs[i]])
		)
	end

	for i = 1:num_piece_types # respect demand
		@constraint(model, edge[i] <= d[i])
	end

	# Not all edges are backward edges, however, we do not have a separate list
	# at the moment. The easiest way is to traverse all edges, if an edge has a
	# non-empty back2indxs, then it has other edges allowing more flow to it and,
	# consequently, it is a backward/circulation edge.
	for i = 1:last_gedge_idx
		@assert i ∉ back2indxs[i]
		if !isempty(back2indxs[i])
			@constraint(model, edge[i] <= sum(edge[back2indxs[i]]))
		end
	end

	# Create fake edges for the pieces.
	for i = 1:num_piece_types
		push!(edges_data, Edge{N, E}(i, zero(N), zero(N), zero(E)))
	end
	# Create the dummy edge for the number of original plates.
	push!(edges_data, Edge{N, E}(
		num_piece_types + 1, zero(N), zero(N), num_piece_types + 2
	))
	# Order the edges by their index.
	sort!(edges_data, by = e -> e.indx)
	#@show num_piece_types
	#@show edges_data[num_piece_types+1:num_piece_types+4]

	BUILT_MODEL, (nodes_data, edges_data)
end # build_model_no_symmbreak

end # module


text:
module Quadratic_Loss
    using LoopVectorization

    function func(output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        loss_matrix = zeros(Float32, size(output_matrix))
        @avxt for i in axes(loss_matrix, 1), j in axes(loss_matrix, 2)
            loss_matrix[i,j] = (output_matrix[i,j]-sample_matrix[i,j])^2
        end
        return loss_matrix
    end

    function prop!(δ::Array{Float32}, output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        @avxt for i in eachindex(i)
            δ[i] = 2*(output_matrix[i]-sample_matrix[i])
        end
    end
end


text:
# This file was generated, do not modify it. # hide
s = schema(X0)
(names=collect(s.names), scitypes=collect(s.scitypes)) |> pretty

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
## Flatten, flatten nested operations that are associative

# eg   (a + b + (c + d)) --> (a + b + c + d)

# FlatT because Flat is already a symbol
#typealias FlatT Union{Mxpr{:Plus},Mxpr{:Times},Mxpr{:And},Mxpr{:Or}, Mxpr{:LCM}, Mxpr{:GCD} }
# We do not rely only on the Flat attribute. We use FlatT in the hope that Julia compiles
# efficient code for each type in the Union.

# Might be faster to interleave the terms because they may need less
# ordering then.  As it is, this is very fast for flattening two long
# sums.  We tried this:
# m1 = Apply(Plus,Table(x^i,[i,1000]));
# m2 = Apply(Plus,Table(x^i,[i,1000])); then
# m3 = mxpr(:Plus,symval(:m1),symval(:m2)) flatten(m3) --> 6e-5s
# canonexpr!(m3) (ie ordering) --> 0.011s running canonexpr! again
# (i.e. it is already sorted and combined) takes 0.00925 s, that is
# only a little bit faster.  So maybe optimizing here is not worth
# anything at the moment.  At Symata cli, m3 = m1 + m2 --> 0.04 s.
# calling setfixed() after canonexpr cuts this time to 0.01s But, this
# cannot be done, in general. Maxima generates the two sums much more
# slowly but adds them much more quickly.

# Flatten one level only
function flatten(mx::T) where T
    needsflat::Bool = false
    for x in margs(mx)
        if isa(x,T)
            needsflat = true
            break
        end
    end
    needsflat == false && return mx
    na = newargs()
    for x in margs(mx)
        if isa(x,T)
            for y in margs(x)
                push!(na,y)
            end
        else
            push!(na,x)
        end
    end
    return mxpr(mhead(mx),na)
end

# Flatten one level only if the argument has the Flat attribute
maybeflatten(mx::FlatT) = flatten(mx)
maybeflatten(x) = x
# Flatten to all levels
function flatten_recursive(mx::T) where T
    needsflat::Bool = false
    for x in margs(mx)
        if isa(x,T)
            needsflat = true
            break
        end
    end
    needsflat == false && return mx
    na = newargs()
    for x in margs(mx)
        if isa(x,T)
            for y in margs(flatten_recursive(x))
                push!(na,y)
            end
        else
            push!(na,x)
        end
    end
    return mxpr(mhead(mx),na)
end

mutable struct FlattenData
    level::Int
    maxlevel::Int
    head::Symbol  # We may want type Any, for heads that are not Symbols
end

# Flatten from level 1 to level n
function flatten_recursive(mx::Mxpr{T}, n::Integer) where T
    d = FlattenData(1,n,T)
    _flatten_recursive(mx,d)
end

# Flatten expressions with head headtype from level 1 to level n
function flatten_recursive(mx::Mxpr, n::Integer, headtype::Symbol)
    d = FlattenData(1,n,headtype)
    _flatten_recursive(mx,d)
end

function _flatten_recursive(mx::Mxpr, d::FlattenData)
    needsflat::Bool = false
    for x in margs(mx)
        if is_Mxpr(x,d.head)
            needsflat = true
            break
        end
    end
    needsflat == false && return mx
    na = newargs()
    for x in margs(mx)
        if is_Mxpr(x,d.head)
            if d.level >= d.maxlevel
                for y in margs(x)
                    push!(na,y)
                end
            else
                d.level += 1
                for y in margs(_flatten_recursive(x,d))
                    push!(na,y)
                end
                d.level -=1
            end
        else
            push!(na,x)
        end
    end
    return mxpr(mhead(mx),na)
end

# TODO: implement the "transpose" case
@mkapprule Flatten :nodefault => true

@sjdoc Flatten """
    Flatten(expr)

remove braces from `Lists` at all levels of `expr`.

    Flatten(expr,n)

flatten only down to level `n`.

    Flatten(expr,n,h)

flatten only expressions with head `h`.
"""

@doap Flatten(x::Mxpr) = flatten_recursive(x)
@doap Flatten(x) = x

@doap function Flatten(x::Mxpr, n::Integer)
    n == 0 && return x
    flatten_recursive(x,n)
end

@doap function Flatten(x::Mxpr, n, headtype::Symbol)
    n == 0 && return x
    n = (n == Infinity ? typemax(Int) : n)
    flatten_recursive(x,n, headtype)
end

@doap function Flatten(x::Mxpr, i::Mxpr{:DirectedInfinity})
    i[1] != 1 && return x
    flatten_recursive(x)
end

@doap function Flatten(x::Mxpr,spec::ListT)
    @warn("Flatten: unimplemented feature")
    return mx
end


text:
"""
Helper functions for distributed Oja's algorithm using MPI

Author: Wayne Wang
Modified: 08/02/2021
"""

function init(d_par::Int, k::Int)
    """
    Initilize eigenvector partition, its gradient partition, and learning rate
    """
    V_par = randn((d_par, k))
    V_par .= Array(qr!(V_par).Q)
    grad_par = similar(V_par)
    α = 1e-5 * ones(k)
    return V_par, grad_par, α
end

function data_par_loader(node_id::Int, input_dir::AbstractString)
    """
    Load a data partition
    """
    X_par = load(string(input_dir, "data_par_", node_id + 1, ".jld2"),
                "data_par")
    return X_par
end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
using Documenter
using ReservoirManagement

makedocs(
    modules = [ReservoirManagement],
    format = :html,
    sitename = "ReservoirManagement.jl",
    authors = "Thibaut Cuvelier.",
    analytics = "UA-67298363-2",
    pages = Any[
        "Home" => "index.md",
        "Scenario generation" => "scenario_generation.md"
    ]
)


text:
using Tabulate
using Base.Test

x1 = ["A","A","B","C","D","D"]
x2 = [true, false, true, false, true, true]
x3 = [3, 5, 4, 2, 2, 3]
x4 = [10.5, 11.7, 13.8, 8.6, 9.2, 6.0]


text:
mutable struct SinglecoefPlot
    thiscoef_loc::Real
    point_est::Real
    confint_lb::Real
    confint_ub::Real
    thiscoef_label::String
    options::SinglecoefplotOption

    function SinglecoefPlot(point_est::Real,
                            confint_lb::Real, confint_ub::Real,
                            thiscoef_label::String,thiscoef_loc::Real,
                            options::SinglecoefplotOption=default_singlecoefplot_options())
        @assert confint_lb <= point_est <= confint_ub "Point estimate not in the interval"
        @assert isfinite(point_est) & isfinite(confint_lb) & isfinite(confint_ub) "Estimation not finite"
        @assert (~ismissing(thiscoef_loc)) || isfinite(thiscoef_loc) "Coef location is not finite"
        new(thiscoef_loc, point_est, confint_lb, confint_ub, thiscoef_label, options)
    end
end

empty_sc(;label::String) = SinglecoefPlot(0, 0, 0, label, 0)

get_line_options(s::SinglecoefPlot) = get_line_options(s.options)
get_dot_options(s::SinglecoefPlot) = get_dot_options(s.options)

function get_dot(singlecoefplot::SinglecoefPlot, vertical::Bool=false)
    loc = (singlecoefplot.thiscoef_loc,singlecoefplot.point_est)
    options = get_dot_options(singlecoefplot)
    if vertical
        Dot(loc, options)
    else
        Dot(reverse(loc), options)
    end
end

function get_line(singlecoefplot::SinglecoefPlot, vertical::Bool=false)
    end1 = (singlecoefplot.thiscoef_loc, singlecoefplot.confint_lb)
    end2 = (singlecoefplot.thiscoef_loc, singlecoefplot.confint_ub)
    options = get_line_options(singlecoefplot)
    if vertical
        LineWithEnds(end1, end2, options)
    else
        LineWithEnds(reverse(end1), reverse(end2), options)
    end
end


text:
using ArrayLayouts, LinearAlgebra, FillArrays, Test
import ArrayLayouts: ApplyBroadcastStyle

@testset "Ldiv" begin
    @testset "Float64 \\ *" begin
        A = randn(5,5)
        b = randn(5)
        M = Ldiv(A,b)

        @test ndims(M) == 1
        @test size(M) == (5,)
        @test similar(M) isa Vector{Float64}
        @test materialize(M) isa Vector{Float64}
        @test all(materialize(M) .=== (A\b) .=== ldiv(A,b))

        @test Base.BroadcastStyle(typeof(Ldiv(A,b))) isa ApplyBroadcastStyle

        @test all(copyto!(similar(b), Ldiv(A,b)) .===
                    (similar(b) .= Ldiv(A,b)) .=== materialize(Ldiv(A,b)) .===
                  (A\b) .=== (b̃ =  copy(b); LAPACK.gesv!(copy(A), b̃); b̃))

        @test copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) ≈ UpperTriangular(A) \ b
        @test all(copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) .===
                    (similar(b) .= Ldiv(UpperTriangular(A),b)) .===
                    BLAS.trsv('U', 'N', 'N', A, b) )

        @test copyto!(similar(b), Ldiv(UpperTriangular(A)' , b)) ≈ UpperTriangular(A)' \ b
        @test all(copyto!(similar(b), Ldiv(UpperTriangular(A)' , b)) .===
                    (similar(b) .= Ldiv(UpperTriangular(A)',b)) .===
                    copyto!(similar(b), Ldiv(transpose(UpperTriangular(A)) , b)) .===
                            (similar(b) .= Ldiv(transpose(UpperTriangular(A)),b)) .===
                    BLAS.trsv('U', 'T', 'N', A, b))
    end

    @testset "ComplexF64 \\ *" begin
        T = ComplexF64
        A = randn(T,5,5)
        b = randn(T,5)
        @test all(copyto!(similar(b), Ldiv(A,b)) .===
                    (similar(b) .= Ldiv(A,b)) .===
                  (A\b) .=== (b̃ =  copy(b); LAPACK.gesv!(copy(A), b̃); b̃))

        @test copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) ≈ UpperTriangular(A) \ b
        @test all(copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) .===
                    (similar(b) .= Ldiv(UpperTriangular(A),b)) .===
                    BLAS.trsv('U', 'N', 'N', A, b) )

        @test copyto!(similar(b), Ldiv(UpperTriangular(A)' , b)) ≈ UpperTriangular(A)' \ b
        @test all(copyto!(similar(b), Ldiv(UpperTriangular(A)' , b)) .===
                    (similar(b) .= Ldiv(UpperTriangular(A)',b)) .===
                    BLAS.trsv('U', 'C', 'N', A, b))

        @test copyto!(similar(b), Ldiv(transpose(UpperTriangular(A)) , b)) ≈ transpose(UpperTriangular(A)) \ b
        @test all(copyto!(similar(b), Ldiv(transpose(UpperTriangular(A)) , b)) .===
                            (similar(b) .= Ldiv(transpose(UpperTriangular(A)),b)) .===
                    BLAS.trsv('U', 'T', 'N', A, b))
    end

    @testset "BigFloat Triangular \\" begin
        for T in (Float64, ComplexF64)
            A = big.(randn(T,5,5))
            b = big.(randn(T,5))
            @test materialize!(Ldiv(UpperTriangular(A),copy(b))) ≈ UpperTriangular(A)\b
            @test materialize!(Ldiv(UnitUpperTriangular(A),copy(b))) ≈ UnitUpperTriangular(A)\b
            @test materialize!(Ldiv(LowerTriangular(A),copy(b))) ≈ LowerTriangular(A)\b
            @test materialize!(Ldiv(UnitLowerTriangular(A),copy(b))) ≈ UnitLowerTriangular(A)\b

            @test materialize!(Ldiv(UpperTriangular(A)',copy(b))) ≈ UpperTriangular(A)'\b
            @test materialize!(Ldiv(UnitUpperTriangular(A)',copy(b))) ≈ UnitUpperTriangular(A)'\b
            @test materialize!(Ldiv(LowerTriangular(A)',copy(b))) ≈ LowerTriangular(A)'\b
            @test materialize!(Ldiv(UnitLowerTriangular(A)',copy(b))) ≈ UnitLowerTriangular(A)'\b

            @test materialize!(Ldiv(transpose(UpperTriangular(A)),copy(b))) ≈ transpose(UpperTriangular(A))\b
            @test materialize!(Ldiv(transpose(UnitUpperTriangular(A)),copy(b))) ≈ transpose(UnitUpperTriangular(A))\b
            @test materialize!(Ldiv(transpose(LowerTriangular(A)),copy(b))) ≈ transpose(LowerTriangular(A))\b
            @test materialize!(Ldiv(transpose(UnitLowerTriangular(A)),copy(b))) ≈ transpose(UnitLowerTriangular(A))\b

            B = big.(randn(T,5,5))
            @test materialize!(Ldiv(UpperTriangular(A),copy(B))) ≈ UpperTriangular(A)\B
        end
    end

    @testset "Triangular \\ matrix" begin
        A = randn(5,5)
        b = randn(5,5)
        M =  Ldiv(UpperTriangular(A), b)
        @test Base.Broadcast.broadcastable(M) === M
        @test UpperTriangular(A) \ b ≈ copyto!(similar(b) , Ldiv(UpperTriangular(A), b)) ≈
            (b .= Ldiv(UpperTriangular(A), b))
    end

    @testset "Int" begin
        A = [1 2 ; 3 4]; b = [5,6];
        @test eltype(Ldiv(A, b)) == Float64
    end

    @testset "Rdiv" begin
        @testset "Float64 \\ *" begin
            A = randn(3,5)
            B = randn(5,5)
            M = Rdiv(A,B)

            @test eltype(M) == Float64
            @test ndims(M) == 2
            @test size(M) == (3,5)
            @test axes(M) == (Base.OneTo(3),Base.OneTo(5))
            @test similar(M) isa Matrix{Float64}
        end
    end

    @testset "Diagonal" begin
        D = Diagonal(randn(5))
        F = Eye(5)
        A = randn(5,5)
        @test copy(Ldiv(D,F)) == ldiv(D,F) == D \ F
        @test ldiv(F,D) == F \ D
        @test ldiv(D,A) ≈ D \ A
        @test ldiv(A,D) == A \ D
        @test ldiv(F,A) ≈ F \ A
        @test ldiv(A,F) == A \ F

        @test copy(Rdiv(D,F)) == rdiv(D,F) == D / F
        @test rdiv(F,D) == F / D
        @test rdiv(A,D) ≈ A / D
        @test rdiv(A,F) == A / F
    end
end


text:
@testset "utils.jl" begin
    rng = MersenneTwister(42)
    m = randn(rng)
    σ = rand(rng)
    d = Normal(m, σ)
    @test AGPL.second_moment(d) ≈ abs2(m) + abs2(σ)
    c = AGPL.second_moment(d)
    @test AGPL.approx_expected_logistic(m, c) ≈ exp(m / 2) / cosh(c / 2) / 2 atol = 1e-5
    for T in [Float64, Float32]
        bigm = T(1000)
        c = bigm + abs(randn(rng, T))
        @test AGPL.approx_expected_logistic(bigm, c) ≈ logistic(bigm)
    end
end


text:
using CSV, DataFrames

base_path = "/Volumes/davidkunz/Documents/bakalarka/mri_head_homography_results_all"
method_folder_names = filter(x -> !occursin(".", x), readdir(base_path))

latexify(df_summary, env=:table, latex=false)

get_mean(name) = round(Float64(results.mean[results.Column1 .== name][1]), sigdigits=3)

function get_method_summary(base_path, method_folder_name)

    method_folder = joinpath(base_path, method_folder_name)
    results_file_path = joinpath(method_folder, "results-summary.csv")
    results = CSV.read(results_file_path) |> DataFrame
    get_mean(name) = Float64(results.mean[results.Column1 .== name][1])
    # rownames = collect(results["Column1"])

    method_name = method_folder_name[findfirst(r"[a-zA-Z]+", method_folder_name)]
    method_name = method_name[1:2] == "Bm" ? method_name[3:end] : method_name
    if method_name == "Sparse"
        open(joinpath(base_path, method_folder_name, "results-summary.txt")) do f
            for i in eachline(f)
                if occursin(r"algorithm", i)
                    method_name = reverse(reverse(i)[findfirst(r"[a-zA-Z_]+", reverse(i))])
                    # @info i method_name reverse(i)
                end
            end
        end
    end
    if method_name == "pflap"
        found = false
        open(joinpath(base_path, method_folder_name, "0", "registration.log")) do f
            for i in eachline(f)
                if occursin(":prefilter => true", i)
                    found = true
                    # @info i method_name reverse(i)
                end
            end
        end
        resized = false
        open(joinpath(base_path, method_folder_name, "results-summary.txt")) do f
            for i in eachline(f)
                if occursin(r"diag_pix", i)
                    resized = true
                    global reduced_size = string(i[findfirst(r"[0-9]+", i)])
                end
            end
        end
        if found
            if resized
                method_name = "\$\\text{PF-LAP}_{\\text{prefilt}}\$\\textsuperscript{*}"
            else
                method_name = "\$\\text{PF-LAP}_{\\text{prefilt}}\$"
            end
        else
            if resized
                method_name = "\$\\text{PF-LAP}\$\\textsuperscript{*}"
            else
                method_name = "\$\\text{PF-LAP}\$"
            end
        end
    elseif method_name == "sparse_pflap"
        resized = false
        open(joinpath(base_path, method_folder_name, "results-summary.txt")) do f
            for i in eachline(f)
                if occursin(r"diag_pix", i)
                    resized = true
                    global reduced_size = string(i[findfirst(r"[0-9]+", i)])
                end
            end
        end
        if resized
            method_name = "\$\\text{Sparse PF-LAP}\\textsuperscript{*}\$"
        else
            method_name = "\$\\text{Sparse PF-LAP}\$"
        end
    elseif method_name == "sparse_pflap_psnr"
        resized = false
        open(joinpath(base_path, method_folder_name, "results-summary.txt")) do f
            for i in eachline(f)
                if occursin(r"diag_pix", i)
                    resized = true
                    global reduced_size = string(i[findfirst(r"[0-9]+", i)])
                end
            end
        end
        if resized
            method_name = "\$\\text{Sparse PF-LAP}_{\\mathrm{PSNR}}\\textsuperscript{*}\$"
        else
            method_name = "\$\\text{Sparse PF-LAP}_{\\mathrm{PSNR}}\$"
        end
    else
        method_name = "\$\\text{$method_name}\$"
    end

    exec_time = "Execution time [minutes]"
    exec_time_secs = "Execution time [seconds]"
    tre_mean = "TRE Mean"
    tre_median = "TRE Median"
    diag_pix = "Image diagonal [pixels]"

    wanted_values = [tre_mean, tre_median, exec_time, diag_pix]

    df_summary = DataFrame(Symbol("Method") => method_name,
                           Symbol("\$\\mathrm{TRE}_{\\text{Mean}}\$") => round(get_mean(tre_mean), sigdigits=3),
                           Symbol("\$\\mathrm{TRE}_{\\text{Med}}\$") => round(get_mean(tre_median), sigdigits=3),
                           Symbol("Time") => round(get_mean(exec_time)*60, sigdigits=3)
                           # Symbol(diag_pix) => get_mean(diag_pix))
                           )
    return df_summary
end

function get_df_all_methods(base_path)
    method_folder_names = filter(x -> !occursin(".", x), readdir(base_path))

    df = get_method_summary(base_path, method_folder_names[1])
    delete!(df, 1:size(df,1))

    for folder_name in method_folder_names
        method_df = get_method_summary(base_path, folder_name)
        # @info eachcol(method_df)
        push!(df, DataFrameRow(method_df, 1))
    end
    return df
end

base_path = "/Volumes/davidkunz/Documents/bakalarka/ox_affine_separate_results/bikes"
base_path = "/Volumes/davidkunz/Documents/bakalarka/ox_affine_separate_results/ubc"
df_out = get_df_all_methods(base_path)

describe(df_out)

base_path = "/Volumes/davidkunz/Documents/bakalarka/mri_head_homography_results_all/"
df_out = get_df_all_methods(base_path)

describe(df_out)

using Latexify

latexify(df_out, env=:table, latex=false)

Dict(eachrow(get_df_summary(base_path, method_folder_names[4])))

## ox affine tables

folders = ["bikes", "leuven", "trees", "ubc"]
base_path = "/Volumes/davidkunz/Documents/bakalarka/ox_affine_separate_results"

base_path = "/Volumes/davidkunz/Documents/bakalarka/"
folders = ["mri_head_homography_results_all"]

for folder in folders
    path = joinpath(base_path, folder)
    df_out = get_df_all_methods(path)
    @info folder
    println(latexify(df_out, env=:table, latex=false))
end


text:
for N in [Float64, Rational{Int}, Float32]

    # π/2 trigonometric rotation
    b = BallInf(N[1, 2], N(1))
    M = N[0 -1 ; 1 0]

    # construction
    lm1 = InverseLinearMap(M, b)
    @test lm1.M == M
    @test lm1.X == b

    # dimension
    @test dim(lm1) == 2

    # Test Support Vector and Support Function
    d = N[1, 1]
    @test σ(d, lm1) == N[3, 0]
    @test ρ(d, lm1) == N(3)
    d = N[-1, 1]
    @test σ(d, lm1) == N[1, 0]
    @test ρ(d, lm1) == N(-1)
    d = N[-1, -1]
    @test σ(d, lm1) == N[1, -2]
    @test ρ(d, lm1) == N(1)
    d = N[1, -1]
    @test σ(d, lm1) == N[3, -2]
    @test ρ(d, lm1) == N(5)

    # scalar multiplication and concretize
    b = BallInf(N[0, 0, 0], N(2))
    ilm1 = InverseLinearMap(N(5), b)
    ilm2 = InverseLinearMap(Matrix(ilm1.M), b)
    @test ilm1 == ilm2
    @test overapproximate(concretize(ilm2.M*ilm2), Hyperrectangle) == overapproximate(b, Hyperrectangle)

    # emptiness check
    @test !isempty(ilm1)
end


text:
# ODEInterface.jl Algorithms

abstract type ODEInterfaceAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
struct dopri5 <: ODEInterfaceAlgorithm end
struct dop853 <: ODEInterfaceAlgorithm end
struct odex <: ODEInterfaceAlgorithm end
struct seulex <: ODEInterfaceAlgorithm end
struct radau <: ODEInterfaceAlgorithm end
struct radau5 <: ODEInterfaceAlgorithm end
struct rodas <: ODEInterfaceAlgorithm end
struct ddeabm <: ODEInterfaceAlgorithm end
struct ddebdf <: ODEInterfaceAlgorithm end


text:
"""
This file contains all the functions that necessary to extract data from GNSS stations that the
Earth Observatory of Singapore (EOS-SG) maintains.

"""

function eosload(reg)

    @info "$(Dates.now()) - Loading information on available GNSS/GPS stations provided by the Earth Observatory of Singapore (EOS-SG)."
    allstn = readdlm(joinpath(@__DIR__,"../data/GNSS-EOS-SG.txt"),',',comments=true);

    @info "$(Dates.now()) - Filtering out for stations in the $(regionfullname(reg)) region."
    lon = allstn[:,2]; lat = allstn[:,3];
    stations = allstn[ispointinregion.(lon,lat,reg),:]; nstns = size(stations,1);

    @info "$(Dates.now()) - There are $(nstns) stations maintained by the Earth Observatory of Singapore (EOS-SG) in the $(regionfullname(reg)) region with Zenith Wet Delay data."

    return stations

end

function eosloadstation(name::AbstractString)

    @info "$(Dates.now()) - Loading information on available GNSS/GPS stations provided by the Earth Observatory of Singapore (EOS-SG)."
    allstn = readdlm(joinpath(@__DIR__,"GNSS-EOS-SG.txt"),',',comments=true);

    @info "$(Dates.now()) - Finding information for station $(name)."
    allnames = allstn[:,1]; stnID = (allnames .== name);
    return allstn[stnID,:]

end

function eosresortsave(zwd::AbstractArray,sig::AbstractArray,
    info::AbstractArray,yrii::Integer,groot::AbstractString)

    ndays = Dates.daysinyear(yrii); nhours = 144;
    zwd = reshape(zwd,nhours,ndays); sig = reshape(sig,nhours,ndays);

    gfol = gnssfol(groot,info[1]); fnc = joinpath(gfol,"$(info[1])-$(yrii).nc");
    if isfile(fnc)
        @info "$(Dates.now()) - Unfinished netCDF file $(fnc) detected.  Deleting."
        rm(fnc);
    end

    ds = Dataset(fnc,"c");
    ds.dim["hour_of_day"] = nhours; ds.dim["day_of_month"] = nlat;

    att_zwd = Dict("units"=>"m","standard_name"=>"zenith_wet_delay",
                   "long_name"=>"Zenith Wet Delay",
                   "scale_factor"=>1/65534,"add_offset"=>0.5,
                   "missing_value"=-32768);
    att_sig = Dict("units"=>"m","standard_name"=>"zenith_wet_delay_error",
                   "long_name"=>"Zenith Wet Delay Uncertainty",
                   "scale_factor"=>0.005/65534,"add_offset"=>0.0025,
                   "missing_value"=-32768);
    att_lon = Dict("units"=>"degrees_east","long_name"=>"longitude");
    att_lat = Dict("units"=>"degrees_north","long_name"=>"latitude");
    att_z   = Dict("units"=>"m","standard_name"=>"height_above_WGS84_surface",
                   "long_name"=>"Orographic Height");

    @info "$(Dates.now()) - Creating GNSS Zenith Wet Delay netCDF file $(fnc) ..."

    vzwd = defVar(ds,"zwd",Int16,("hour_of_day","day_of_month"),attrib=att_zwd);
    vsig = defVar(ds,"sig",Int16,("hour_of_day","day_of_month"),attrib=att_sig);
    vzwd.var[:] = zwd; vsig.var[:] = sig;

    defVar(ds,"longitude",info[2],attrib=att_lon);
    defVar(ds,"latitude",info[3],attrib=att_lat);
    defVar(ds,"z",info[4],attrib=att_z);

    close(ds);

    @info "$(Dates.now()) - Zenith Wet Delay data for $(info[1]) has been saved into file $(fnc)."

end

function eosextractyear(gnssdata::AbstractArray,info::AbstractArray,
    yrArray::AbstractArray,yrii::Integer)

    yrdates = convert(Array,DateTime(yrii,1,1):Minute(10):DateTime(yrii+1,1,1)); pop!(yrdates);

    @debug "$(Dates.now()) - Extracting GNSS Zenith Wet Delay data for Year $(yrii) ..."
    gnssdata = gnssdata[findall(isequal(yrii),yrArray),:];
    zwd = zeros(size(yrdates)); sig = zeros(size(yrdates)); ei = 0;

    for ti = 1 : size(yrdates,1)
        try jj = findfirst(isequal(yrdates[ti]),gnssdata[:,2])
               zwd[ti] = gnssdata[jj,3]; sig[ti] = gnssdata[jj,4]
        catch; zwd[ti] = NaN; sig[ti] = NaN; ei = ei + 1;
        end
    end

    @debug "$(Dates.now()) - There are $(ei) entries without any valid data for Year $(yrii)."
    return zwd,sig

end

function eosresort(stations::AbstractArray,groot::Dict)

    for ii = 1 : size(stations,1)

        fraw  = joinpath(groot["raw"],"$(stations[ii,1]).tdpzwd");
        info  = stations[ii,:];

        @info "$(Dates.now()) - Available data from $(stations[ii,1]) GNSS station will be extracted from file $(fraw)."
        gdata = readdlm(fraw,Any,comments=true);

        @info "$(Dates.now()) - Converting J2000 seconds data into Julia DateTime format ..."
        gdata[:,2] .= DateTime(2000,1,1,12,0,0) + Second.(gdata[:,2]);
        yr = Dates.year.(gdata[:,2]); yrbeg = minimum(yr); yrend = maximum(yr);

        @info "$(Dates.now()) - Extracting Zenith Wet Delay data from the GNSS station $(stations[ii,1]) ..."
        for yrii = yrbeg : yrend
            zwd,sig = eosextractyear(gdata,info,yr,yrii);
            eosresortsave(zwd,sig,info,yrii,groot["data"]);
        end
        @info "$(Dates.now()) - Zenith Wet Delay data from the GNSS station $(stations[ii,1]) has been extracted and saved into yearly NetCDF files."

    end

end