text:
using SolveDSGE

filename = "model1b.txt"
path = joinpath(@__DIR__,filename)
process_model(path)
dsge = retrieve_processed_model(path)

x = [0.05, 21.4, 1.97, 2.8]

tol = 1e-8
maxiters = 1000
ss_obj = compute_steady_state(dsge,x,tol,maxiters)
ss = ss_obj.zero

N   = PerturbationScheme(ss,1.0,"first")
NN  = PerturbationScheme(ss,1.0,"second")
NNN = PerturbationScheme(ss,1.0,"third")

soln_fo = solve_model(dsge,N)
soln_so = solve_model(dsge,NN)
soln_to = solve_model(dsge,NNN)

P = ChebyshevSchemeDet(ss,chebyshev_nodes,[21,21],3,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
PP = ChebyshevSchemeDet(ss,chebyshev_nodes,[21,21],4,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
PPP = ChebyshevSchemeDet(ss,chebyshev_nodes,[71,71],6,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)

soln_nla = solve_model(dsge,P)
soln_nlb = solve_model(dsge,soln_fo,PP)
soln_nlc = solve_model(dsge,soln_so,PPP,2)
soln_nld = solve_model(dsge,soln_to,P)
soln_nle = solve_model(dsge,soln_nld,PP)
soln_nlf = solve_model(dsge,soln_nle,PPP,2)

L = SmolyakSchemeDet(ss,chebyshev_gauss_lobatto,4,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
LL = SmolyakSchemeDet(ss,clenshaw_curtis_equidistant,4,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
LLL = SmolyakSchemeDet(ss,chebyshev_gauss_lobatto,5,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)

soln_nlg = solve_model(dsge,L)
soln_nlh = solve_model(dsge,soln_to,LL)
soln_nli = solve_model(dsge,soln_nlh,LLL,2)
soln_nlj = solve_model(dsge,soln_nlf,LL)

M = PiecewiseLinearSchemeDet(ss,[21,21],[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
MM = PiecewiseLinearSchemeDet(ss,[31,31],[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)

soln_nlk = solve_model(dsge,M)
soln_nll = solve_model(dsge,soln_to,MM)
soln_nlm = solve_model(dsge,soln_nlf,MM,2)
soln_nln = solve_model(dsge,soln_nlj,MM)
soln_nlo = solve_model(dsge,soln_nln,MM)

simulated_data1 = simulate(soln_nlf,0.95*ss[1:2],100000)


text:
"""
Returns the Finite Impulse Response of a Raised Cosine (RC) filter.	The filter is defined by its span (evaluated in number of symbol N), its Roll-Off factor and its oversampling factor. The span corresponds to the number of symbol affected by filter before and after the center point.\n
Output is a Vector{Float64} array of size L= 2KN+1 \n
SRRC definition is based on [1]	 \n
[1]	 3GPP TS 25.104 V6.8.0 (2004-12). http://www.3gpp.org/ftp/Specs/archive/25_series/25.104/25104-680.zip \n
 Syntax \n
		h	= raisedCosine(N,beta,ovS)
Input parameters \n
- N	  : Symbol span (Int16)
- beta  : Roll-off factor (Float64)
- ovS	  : Oversampling rate (Int16)
"""
function raisedCosine(N,beta,ovS)
	# --- Final size of filter
	nbTaps	= 2 * N * ovS + 1;
	# --- Init output
	h			= zeros(Float64,nbTaps);
	counter		= 0;
	# --- Iterative SRRC definition
	for k = -N*ovS : 1 : N*ovS
		counter		 = counter + 1;
		if k == 0
			## First singular point at t=0
			h[counter]  = 1;#(1-beta) + 3*beta/pi;
		elseif abs(k) == ovS / (2*beta);
			## Second possible singular point
			h[counter]  = pi/4*sin(pi/(2beta))/(pi/(2beta));
		else
			## Classic SRRC formulation (see [1])
			h[counter]	= sin(pi*k/ovS)/(pi*k/ovS) * cos(pi*beta*k/ovS)/(1- 4beta^2*(k/ovS)^2);
		end
	end
	return h
end

"""
Returns the Finite Impulse Response of a Square Root Raised Cosine (SRRC) filter. \n
The filter is defined by its span (evaluated in number of symbol N), its Roll-Off factor and its oversampling factor. The span corresponds to the number of symbol affected by filter before and after the center point.\n
Output is a Vector{Float64} array of size L= 2KN+1\n
SRRC definition is based on [1]\n
[1]	 3GPP TS 25.104 V6.8.0 (2004-12). http://www.3gpp.org/ftp/ Specs/archive/25_series/25.104/25104-680.zip\n
Syntax\n
h	= sqrtRaisedCosine(N,beta,ovS) \n
Input parameters \n
- N	  : Symbol span (Int16)
- beta  : Roll-off factor (Float64)
- ovS	  : Oversampling rate (Int16)
"""
function sqrtRaisedCosine(N,beta,ovS)
	# --- Final size of filter
	nbTaps	= 2 * N * ovS + 1;
	# --- Init output
	h			= zeros(Float64,nbTaps);
	counter		= 0;
	# --- Iterative SRRC definition
	for k = -N*ovS : 1 : N*ovS
		counter		 = counter + 1;
		if k == 0
			## First singular point at t=0
			h[counter]  = (1-beta) + 4*beta/pi;
		elseif abs(k) == ovS / (4*beta);
			## Second possible singular point
			h[counter]  = beta/sqrt(2)*( (1+2/pi)sin(pi/(4beta))+(1-2/pi)cos(pi/(4beta)));
		else
			## Classic SRRC formulation (see [1])
			h[counter]  = ( sin(pi*k/ovS*(1-beta)) + 4beta*k/ovS*cos(pi*k/ovS*(1+beta))) / (pi*k/ovS*(1- (4beta*k/ovS)^2) );
		end
	end
    # --- Max @ h[0]
    h   = h ./ maximum(h)
	return h
end

ext: jl


text:
d=2;
A=Array{Int64,3}(d,d,d); B=Array{Int64,3}(d,d,d); C=zeros(Int64,d,d)
for i in 1:d,j in 1:d,k in 1:d
    A[i,j,k] = 2i+j^2-k^3
    B[i,j,k] = i+2j+4k
end

tic()
for i in 1:size(A,1)
    for j in 1:size(A,3)
        for k in 1:size(A,1)
            for n in 1:size(B,3)
                C[j,n] += A[i,j,k]*B[k,i,n]
            end
        end
    end
end
toc()
print("Matrix C is ", C)

A_permuted = permutedims(A,[2,1,3]); B_permuted = permutedims(B,[2,1,3]);
A_reshaped = reshape(A_permuted,2,4); B_reshaped = reshape(B_permuted,4,2);
tic()
C_reshaped = A_reshaped*B_reshaped;
toc()
print("Matrix C is ", C_reshaped)

A=rand(Int64,100,100,100); B=rand(Int64,100,100,100);
C=zeros(Int64,100,100);

tic()
for i in 1:size(A,1)
    for j in 1:size(A,3)
        for k in 1:size(A,1)
            for n in 1:size(B,3)
                C[j,n] += A[i,j,k]*B[k,i,n]
            end
        end
    end
end
toc()

A_reshaped=reshape(permutedims(A,[2,1,3]),100,10000);
B_reshaped=reshape(permutedims(B,[2,1,3]),10000,100);
tic()
C_reshaped = A_reshaped*B_reshaped;
toc()

X=Complex{Float64}[0.0 1.0;1.0 0.0];
Y=Complex{Float64}[0.0im -1.0im;1.0im 0.0im];
Z=Complex{Float64}[1.0 0.0;0.0 -1.0];
h=kron(X,X)+kron(Y,Y)+kron(Z,Z);
N=4; C=2; d=2;
H=zeros(Complex{Float64},d^N,d^N);

for l in 0:N-C
    H+=kron(kron(eye(d^l),h),eye(d^(N-C-l)))
end
E,V = eig(H)
print("The first eigenvalues are: ", E[1:11])

function diagonalization_model(h::Matrix{T},N::Int64;C::Int64=2,
        d::Int64=2,method_sparse::Bool=true) where {T}
    H=zeros(T,d^N,d^N)
    for l in 0:N-C
        H+=kron(kron(eye(d^l),h),eye(d^(N-C-l),d^(N-C-l)))
    end
    if method_sparse
        E,V = eigs(H;nev=10,which=:SR)
    else
        E,V = eig(H)
    end
    return E
end

@time diagonalization_model(h,12)

@time diagonalization_model(h,12,method_sparse=false)

function permR2(i::Int64,P::Int64)
    return mod.(collect(1:P)+i-2,P)+1
end

function permL2(i::Int64,P::Int64)
    return mod.(collect(1:P)-i,P)+1
end

function H_wavefunction(h_tilde::Matrix{T},phi::Array{T,P},N::Int64;
        OBC::Int64=1,d::Int64=2,C::Int64=2) where {T,P}
    H_phi = zeros(T,fill(d,C+1)...,d^(N-C-1))
    for i in 1:N-OBC
        if N == 1
            phi = permutedims(phi,permR2(1,P))
            H_phi = permutedims(H_phi,permR2(1,P))
            phi_tmp = h_tilde*reshape(phi,d^C,d^(N-C))
            phi_tmp = reshape(phi_tmp,fill(d,C+1)...,d^(N-C-1))
            H_phi = reshape(H_phi,fill(d,C+1)...,d^(N-C-1))
            H_phi += phi_tmp
        else
            phi = permutedims(phi,permR2(2,P))
            H_phi = permutedims(H_phi,permR2(2,P))
            phi_tmp = h_tilde*reshape(phi,d^C,d^(N-C))
            phi_tmp = reshape(phi_tmp,fill(d,C+1)...,d^(N-C-1))
            H_phi = reshape(H_phi,fill(d,C+1)...,d^(N-C-1))
            phi = reshape(phi,fill(d,C+1)...,d^(N-C-1))
            H_phi += phi_tmp
        end
    end
    if OBC == 1
        return(reshape(permutedims(H_phi,permR2(2,P)),d^N))
    elseif OBC == 0
        return(reshape(H_phi,d^N))
    end
end

function power_method(h::Matrix{T},N::Int64;OBC::Int64=1,d::Int64=2,
        C::Int64=2,tol::Float64=1e-4) where {T}
    #phi = rand(T,fill(d,C+1)...,d^(N-C-1))
    #phi = collect(Complex{Float64},1:2^N)
    phi = rand(T,d^N)
    phi = phi/sqrt(dot(conj(phi),phi))
    phi = reshape(phi,fill(d,C+1)...,d^(N-C-1))
    list_E_tmp = []
    E,V = eig(h)
    alpha=max(E...)
    h_tilde = h - alpha*kron(eye(d),eye(d))
    i = 1
    bool = true
    while bool
        H_phi = H_wavefunction(h_tilde,phi,N,OBC=OBC)
        #H_phi = H_phi/dot(conj(H_phi),H_phi)
        phi = reshape(phi,d^N)
        #E = dot(conj(phi),H_phi)/dot(conj(phi),phi) + α*(N-OBC)
        E = dot(conj(phi),H_phi) + alpha*(N-OBC)
        #println("E = ", E)
        push!(list_E_tmp,E)
        if i>2 && abs(list_E_tmp[i]-list_E_tmp[i-1])<tol
           bool = false
           println("Desired precision reached!")
        end
        i+=1
        phi = reshape(phi,fill(d,C+1)...,d^(N-C-1))
        phi = reshape(H_phi/sqrt(dot(conj(H_phi),H_phi)),fill(d,C+1)...,d^(N-C-1))
    end
    println("E = ",E)
    println("Number of iterations = ", i)
end

@time power_method(h,4,tol=1e-9)

@time power_method(h,12,tol=1e-9)

@time power_method(h,16,OBC=0,tol=1e-9)

@time power_method(h,16,tol=1e-9)

@time power_method(h,25)

@time power_method(h,25,OBC=0)

max_stars_repo_name: oliviersimard/SuperStiff

max_issues_repo_name: oliviersimard/SuperStiff

max_forks_repo_name: oliviersimard/SuperStiff


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    AbstractMCMCTuningStrategy

Abstract super-type for MCMC tuning strategies.
"""
abstract type AbstractMCMCTuningStrategy end
export AbstractMCMCTuningStrategy

abstract type AbstractMCMCTuner end

function isvalid end
function isviable end

function tuning_init! end

function mcmc_tune_burnin! end

function mcmc_init end

"""
    @with_kw struct MCMCBurninStrategy

Defines the MCMC burn-in strategy, specifically the number and length of
MCMC tuning/burn-in cycles.

Fields:

* `max_nsamples_per_cycle`: Maximum number of MCMC samples to generate per
  cycle, defaults to `1000`. Definition of a sample depends on MCMC algorithm.

* `max_nsteps_per_cycle`: Maximum number of MCMC steps per cycle, defaults
  to `10000`. Definition of a step depends on MCMC algorithm.

* `max_time_per_cycle`: Maximum wall-clock time to spend per cycle, in
  seconds. Defaults to `Inf`.

* `max_ncycles`: Maximum number of cycles.
"""
@with_kw struct MCMCBurninStrategy
    max_nsamples_per_cycle::Int64 = 1000
    max_nsteps_per_cycle::Int64 = 10000
    max_time_per_cycle::Float64 = Inf
    max_ncycles::Int = 30
end

export MCMCBurninStrategy

function MCMCBurninStrategy(algorithm::MCMCAlgorithm, nsamples::Integer, max_nsteps::Integer, tuner_config::AbstractMCMCTuningStrategy)
    max_nsamples_per_cycle = max(div(nsamples, 10), 10)
    max_nsteps_per_cycle = max(div(max_nsteps, 10), 10)
    MCMCBurninStrategy(
        max_nsamples_per_cycle = max_nsamples_per_cycle,
        max_nsteps_per_cycle = max_nsteps_per_cycle
    )
end

function mcmc_tune_burnin!(
    callbacks,
    tuners::AbstractVector{<:AbstractMCMCTuner},
    chains::AbstractVector{<:MCMCIterator},
    convergence_test::MCMCConvergenceTest,
    burnin_strategy::MCMCBurninStrategy;
    strict_mode::Bool = false
)
    @info "Begin tuning of $(length(tuners)) MCMC chain(s)."

    nchains = length(chains)

    user_callbacks = mcmc_callback_vector(callbacks, eachindex(chains))

    cycles = zero(Int)
    successful = false
    while !successful && cycles < burnin_strategy.max_ncycles
        stats = [x.stats for x in tuners] # ToDo: Find more generic abstraction

        # Clear all stats before tuning cycle
        empty!.(stats)

        cycles += 1
        run_tuning_cycle!(
            user_callbacks, tuners, chains;
            max_nsamples = burnin_strategy.max_nsamples_per_cycle,
            max_nsteps = burnin_strategy.max_nsteps_per_cycle,
            max_time = burnin_strategy.max_time_per_cycle
        )

        ct_result = check_convergence!(convergence_test, chains, stats)

        ntuned = count(c -> c.info.tuned, chains)
        nconverged = count(c -> c.info.converged, chains)
        successful = (ntuned == nconverged == nchains)

        for i in eachindex(user_callbacks, tuners)
            user_callbacks[i](1, tuners[i])
        end

        @info "MCMC Tuning cycle $cycles finished, $nchains chains, $ntuned tuned, $nconverged converged."
    end

    if successful
        @info "MCMC tuning of $nchains chains successful after $cycles cycle(s)."
    else
        msg = "MCMC tuning of $nchains chains aborted after $cycles cycle(s)."
        if strict_mode
            @error msg
        else
            @warn msg
        end
    end

    successful
end

struct NoOpTunerConfig <: BAT.AbstractMCMCTuningStrategy end

(config::NoOpTunerConfig)(chain::MCMCIterator; kwargs...) = NoOpTuner()

struct NoOpTuner{C<:MCMCIterator} <: AbstractMCMCTuner end

isvalid(chain::MCMCIterator) = current_sample(chain).logdensity > -Inf

isviable(tuner::NoOpTuner, chain::MCMCIterator) = true

function tuning_init!(tuner::NoOpTuner, chain::MCMCIterator)
    nothing
end

function mcmc_tune_burnin!(
    callbacks,
    tuners::AbstractVector{<:NoOpTuner},
    chains::AbstractVector{<:MCMCIterator},
    convergence_test::MCMCConvergenceTest,
    burnin_strategy::MCMCBurninStrategy;
    kwargs...
)
    @debug "Tune/Burn-In with NoOpTuner doing nothing."
end

"""
    @with_kw struct MCMCInitStrategy

Defines the MCMC chain initialization strategy.

Fields:

* `init_tries_per_chain`: Interval that specifies the minimum and maximum
  number of tries per MCMC chain to find a suitable starting position. Many
  candidate chains will be created and run for a short time. The chains with
  the best performance will be selected for tuning/burn-in and MCMC sampling
  run. Defaults to `IntervalSets.ClosedInterval(8, 128)`.

* `max_nsamples_init`: Maximum number of MCMC samples for each candidate
  chain. Defaults to 25. Definition of a sample depends on sampling algorithm.

* `max_nsteps_init`: Maximum number of MCMC steps for each candidate chain.
  Defaults to 250. Definition of a step depends on sampling algorithm.

* `max_time_init::Int`: Maximum wall-clock time to spend per candidate chain,
  in seconds. Defaults to `Inf`.
"""
@with_kw struct MCMCInitStrategy
    init_tries_per_chain::ClosedInterval{Int64} = ClosedInterval(8, 128)
    max_nsamples_init::Int64 = 25
    max_nsteps_init::Int64 = 250
    max_time_init::Float64 = Inf
end

export MCMCInitStrategy

MCMCInitStrategy(tuner_config::AbstractMCMCTuningStrategy) =
    MCMCInitStrategy()

_gen_chains(ids::AbstractRange{<:Integer},
    chainspec::MCMCSpec,
    ) = [chainspec(id) for id in ids]

function mcmc_init(
    chainspec::MCMCSpec,
    nchains::Int,
    tuner_config::AbstractMCMCTuningStrategy = AbstractMCMCTuningStrategy(algorithm(chainspec)),
    init_strategy::MCMCInitStrategy = MCMCInitStrategy(tuner_config)
)
    @info "Trying to generate $nchains viable MCMC chain(s)."

    min_nviable::Int = minimum(init_strategy.init_tries_per_chain) * nchains
    max_ncandidates::Int = maximum(init_strategy.init_tries_per_chain) * nchains

    ncandidates::Int = 0

    dummy_chain = chainspec(zero(Int64))
    dummy_tuner = tuner_config(dummy_chain)

    chains = similar([dummy_chain], 0)
    tuners = similar([dummy_tuner], 0)
    cycle::Int = 1

    while length(tuners) < min_nviable && ncandidates < max_ncandidates
        n = min(min_nviable, max_ncandidates - ncandidates)
        @debug "Generating $n $(cycle > 1 ? "additional " : "")MCMC chain(s)."

        new_chains = _gen_chains(ncandidates .+ (one(Int64):n), chainspec)

        filter!(isvalid, new_chains)

        new_tuners = tuner_config.(new_chains)
        tuning_init!.(new_tuners, new_chains)
        ncandidates += n

        @debug "Testing $(length(new_tuners)) MCMC chain(s)."

        # ToDo: Use mcmc_iterate! instead of run_tuning_iterations! ?
        run_tuning_iterations!(
            (), new_tuners, new_chains;
            max_nsamples = max(5, div(init_strategy.max_nsamples_init, 5)),
            max_nsteps =  max(50, div(init_strategy.max_nsteps_init, 5)),
            max_time = init_strategy.max_time_init / 5
        )

        viable_idxs = findall(isviable.(new_tuners, new_chains))
        viable_tuners = new_tuners[viable_idxs]
        viable_chains = new_chains[viable_idxs]

        @debug "Found $(length(viable_idxs)) viable MCMC chain(s)."

        if !isempty(viable_tuners)
            # ToDo: Use mcmc_iterate! instead of run_tuning_iterations! ?
            run_tuning_iterations!(
                (), viable_tuners, viable_chains;
                max_nsamples = init_strategy.max_nsamples_init,
                max_nsteps = init_strategy.max_nsteps_init,
                max_time = init_strategy.max_time_init
            )

            nsamples_thresh = floor(Int, 0.8 * median([nsamples(chain) for chain in viable_chains]))
            good_idxs = findall(chain -> nsamples(chain) >= nsamples_thresh, viable_chains)
            @debug "Found $(length(viable_tuners)) MCMC chain(s) with at least $(nsamples_thresh) samples."

            append!(chains, view(viable_chains, good_idxs))
            append!(tuners, view(viable_tuners, good_idxs))
        end

        cycle += 1
    end

    length(tuners) < min_nviable && error("Failed to generate $min_nviable viable MCMC chains")

    m = nchains
    tidxs = LinearIndices(tuners)
    n = length(tidxs)

    mode_1 = tuners[1].stats.mode
    modes = Array{eltype(mode_1)}(undef, length(mode_1), n)
    for i in tidxs
        modes[:,i] = tuners[i].stats.mode
    end

    final_chains = similar(chains, zero(Int))
    final_tuners = similar(tuners, zero(Int))

    if 2 <= m < size(modes, 2)
        clusters = kmeans(modes, m)
        clusters.converged || error("k-means clustering of MCMC chains did not converge")

        mincosts = fill(Inf, m)
        chain_sel_idxs = fill(0, m)

        for i in tidxs
            j = clusters.assignments[i]
            if clusters.costs[i] < mincosts[j]
                mincosts[j] = clusters.costs[i]
                chain_sel_idxs[j] = i
            end
        end

        @assert all(j -> j in tidxs, chain_sel_idxs)

        for i in sort(chain_sel_idxs)
            push!(final_chains, chains[i])
            push!(final_tuners, tuners[i])
        end
    else
        @assert length(chains) == nchains
        resize!(final_chains, nchains)
        copyto!(final_chains, chains)

        @assert length(tuners) == nchains
        resize!(final_tuners, nchains)
        copyto!(final_tuners, tuners)
    end

    @info "Selected $(length(final_tuners)) MCMC chain(s)."

    (chains = final_chains, tuners = final_tuners)
end


text:
module Futbolista

using DataFrames
using JSON: Parser, parsefile
using HTTP
using Statistics
using Plots

include("DataLoaders.jl")
include("Plotter.jl")

using .Plotter
using .DataLoaders

end ##end of module


text:
# tests of basic utilities for working with HTML

import Gumbo: HTMLNode, NullNode

# convenience constructor works
@test HTMLElement(:body) == HTMLElement{:body}(HTMLNode[],
                                               NullNode(),
                                               Dict{String,String}())

# accessing tags works
@test HTMLElement(:body) |> tag == :body


text:
# This file was generated, do not modify it. # hide
using MLJ
using StableRNGs
import DataFrames
MLJ.color_off() # hide
@load RidgeRegressor pkg=MultivariateStats

rng = StableRNG(6616) # for reproducibility
x1 = rand(rng, 300)
x2 = rand(rng, 300)
x3 = rand(rng, 300)
y = exp.(x1 - x2 -2x3 + 0.1*rand(rng, 300))
X = DataFrames.DataFrame(x1=x1, x2=x2, x3=x3)

test, train = partition(eachindex(y), 0.8);


text:
"""
    Orthogonal(U, S; validate_fields)

An AbstractMatrix `H` that takes the form `H = U * sqrt(S)` with `U` a matrix
with orthonormal columns and `S` a diagonal matrix with positive entries.

# Arguments:
- U: a `p x m` matrix with mutually orthonormal columns.
- S: an `m x m` `Diagonal` matrix with positive entries.
"""
struct Orthogonal{T<:Real, TU<:AbstractMatrix{T}, TS<:Diagonal{T}} <: AbstractMatrix{T}
    U::TU
    S::TS

    function Orthogonal(U::AbstractMatrix{T}, S::Diagonal{T}; validate_fields=true) where T<:Real
        validate_fields && _validate(U)
        return new{T, typeof(U), typeof(S)}(U, S)
    end
end

function _validate(U::AbstractMatrix)
    isapprox(U' * U, I) || throw(ArgumentError("`U` is not an orthogonal matrix"))
end

@non_differentiable _validate(::AbstractMatrix)

function Base.getindex(h::Orthogonal, args...)
    H = h.U * sqrt(h.S)
    return getindex(H, args...)
end

function Base.size(h::Orthogonal)
    return size(h.U)
end


text:
using KentFinancialEconomics
using Documenter

DocMeta.setdocmeta!(KentFinancialEconomics, :DocTestSetup, :(using KentFinancialEconomics); recursive=true)

makedocs(;
    modules=[KentFinancialEconomics],
    authors="alfredjmduncan <alfie.j.duncan@gmail.com> and contributors",
    repo="https://github.com/alfredjmduncan/KentFinancialEconomics.jl/blob/{commit}{path}#{line}",
    sitename="KentFinancialEconomics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://alfredjmduncan.github.io/KentFinancialEconomics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/alfredjmduncan/KentFinancialEconomics.jl",
)

max_stars_repo_name: alfredjmduncan/KentFinancialEconomics

max_issues_repo_name: alfredjmduncan/KentFinancialEconomics

max_forks_repo_name: alfredjmduncan/KentFinancialEconomics


text:
"""
at this point we already have res list filled with all fp and fn's in some cases it is enough in others we want futher actions to be done
  - for example in case of lack of isometry in voxels we can now apply some corrections for it
"""
module ResultProcess

"""
it will scan through result list - ignorea all that are zeroces in res indicies, add only the iteration numbers and in the end we will just dividy it by (fp+fn)
entriesPerBlock - amount of entries per block of threads
totalLength - total length of result list
iterLoopResList - how many times given block need to iterate to aalyze all entries given to it
globalSum - global variable accessed atomically holding sum of iter numbs of results
"""
function getAverage(resList,resListIndicies,entriesPerBlock,totalLength,iterLoopResList ,globalSum )
 locSum = UInt32(0)
 shmemSum = @cuStaticSharedMem(UInt32, (33,1))
 offsetIter= UInt8(1)

  @iterateLinearlyForResList(iterLoopResList,entriesPerBlock,totalLength,begin
      if(resListIndicies[i]>0)# we do not want 0 entries or duplicated ones
        locSum+=resList[i,6]#adding iteration number
      end
  end )

  #time for reduction
  @redWitAct(offsetIter,shmemSum, locSum,+   )
  sync_threads()
  @addAtomic(shmemSum, globalSum)
end #getAverage

"""
iterate linearly result list
entriesPerBlock - amount of entries per block of threads
totalLength - total length of result list
"""
macro iterateLinearlyForResList(iterLoop,entriesPerBlock,totalLength, ex)
  return  esc(quote
  i = UInt32(0)
  @unroll for j in 0:($iterLoop)
    offset = ((blockIdxX()-1) *$entriesPerBlock)
    i= threadIdxX()+(threadIdxY()-1)*blockDimX()+ j* blockDimX()*blockDimY()
    if((i+offset)<=$totalLength && i<$entriesPerBlock)
      i+=offset
      $ex
    end
  end
   end)
  end

"""
We will iterate over result list take all entries that are non zeo in res  indicies list now in a while loop we will lok for the source voxel that was responsible for the  result
so we can start fom naive approach where each lane will be responsible for finding source voxel that is the closest voxel from other mask
In order to achieve it we will go in the direction indicated by the dir value in the distance indicated by the dir in the result matrix  and the number of voxels indicated by the iteration number
It is probable that we will not find there he source voxel so we will need to iteratively increase number of voxels that go in orthogonal direction so for example
  when direction is down and iteration number 10 we will ook 10 voxels blow in straight line at first
  then if wew ould not find we will look into voxel 9 down and in 1 left,right,anterior,posterior - so over orthogonal axes
  then we will  try going 8 pixels down and 1 pixel in 4 directions and futher pixel also in 4 directions  so 16 permutations
  basically we will want to scan the area of the sphere with radius of iteration number and getting circle (or sth more like a square/rhombus as we analze 1 norm ...) crossections of it  starting from the point in the end of the axis in given
  direction and progress up circle crossesction by circle crossesction
so the amount of voxels that we have to use at each iteration will increase in each step hence w have no free voxels to move at first step - then we have 4 possible position then I suppose 16...
so we can do it like that as we are fro example on step 3 and we are circling around with sth like circle around with radius 3 we can first analyze one axis for example left right so we go
3 voxls laeft and right we have a part of line now we get shorter line (by 1 on each end) anteriorly and posteriorly to the first axis then next row even shorter anteriorly and posteriorly ...
until we will reduce row length to 1 - this should lead to check the voxels at the end of anterior posterior ends of anterior posterior axis

We can xperiment with getting corrections one by single line at the begining with some maximum number of ccles allowed - so if single lane wuld not be able to find result we would cease looking for it
and then scan one more time whole warps  for scanning those points that were not yet resolved

additionally if we want to create intresting visualization we can add information about corrected distance to all points between source and target voxel..

  entriesPerBlock - amount of entries per block of threads
  totalLength - total length of result list
  iterLoopResList - how many times given block need to iterate to aalyze all entries given to it
  globalSum - global variable accessed atomically holding sum of iter numbs of results - so we will accumulate corrected values - to later return average
  maxSingleThrIterNumb - maximum iteration number that should be analyzed using just a single lane; all points bigger than that will be analyzed used whole warps or maybe blocks
  """
function applyCorrection(resList,resListIndicies,entriesPerBlock,totalLength,iterLoopResList,globalSum,maxSingleThrIterNumb, referenceArrs )
  locSum = UInt32(0)
  shmemSum = @cuStaticSharedMem(UInt32, (33,1))
  offsetIter= Int16(1)
  locIterNumb=Int16(0)
  #indicates how coordinates should change relatively to the covered point from result list
  xChange= Int16(0)
  yChange= Int16(0)
  zChange= Int16(0)

@iterateLinearlyForResList(iterLoopResList,entriesPerBlock,totalLength,begin
  if(resListIndicies[i]>0)# we do not want to analyze 0 entries or duplicated ones

  end
end )

end

### looking for the true entry in given array  we iteratively scan area  where the original point from which dilatation started and led to this result that we are currently analyzing

"""
using single thread we iteratively scan area  where the original point from which dilatation started and led to this result that we are currently analyzing
referenceArrs= (reducedGoldB,reducedSegmB)
                top 6
                bottom 5
                left 2
                right 1
                anterior 3
                posterior 4
"""
macro singleThreadScan(resList, i, referenceArrs )

  return  esc(quote

  offsetIter=$resList[$i,6] #storing iteration number
  locIterNumb=0 #will get incemented every time we move centrally in main axis
  while offsetIter>0
    #first we need to  modify x,y or z depending on the direction given in dir variable of result list
    dir = $resList[$i,5]
    xChange-= (dir==2)*offsetIter
    xChange+= (dir==1)*offsetIter

    yChange-= (dir==4)*offsetIter
    yChange+= (dir==3)*offsetIter

    zChange-= (dir==6)*offsetIter
    zChange+= (dir==5)*offsetIter

    #now we have established position over main axis - what is left is to scan by modifying orthogonal axes

    @unroll for orthoAxAmove in -locIterNumb:locIterNumb#so we need to iterate over one of the
      @unroll for orthoAxBmove in [-locIterNumb+ CUDA.abs(orthoAxAmove),locIterNumb-CUDA.abs(orthoAxAmove)]
        #isXMain,isYMain - true if x or in second case y is main axis - we define main axis on the basis of dir variable
        isXMain = (dir==2 || dir==1)
        isYMain = (dir==3 || dir==4)

        #(dir!=2 && dir!=1) so x axis is not main axis
        xChange+= (!isXMain)*orthoAxAmove
        #(dir==2 || dir==1)*orthoAxBmove - so if x is main axis - it means that y is not  and it will be axis A
        yChange+=(!isXMain && !isYMain)*orthoAxBmove+ (isXMain)*orthoAxAmove
        #zChange can never affect axis A  will affect axis B if it is not main axis - so when either x or y is main axis
        zChange+=(!isXMain && !isYMain)*orthoAxBmove

        #here if we will find that there is true in this position we will assume  that those coordinates are coordinates of
        if(referenceArrs[2-$resList[$i,4]][ $resList[$i,1]+xChange,$resList[$i,2]+yChange,$resList[$i,3]+zChange])
            #now on the basis of new coordinates we need to calculate correction for non isometric voxels
            corrected =getCorrectedDistance(yDimSize,zDimSize, $resList[$i,1]+xChange,$resList[$i,2]+yChange,$resList[$i,3]+zChange)
            locSum+=corrected
            $resList[$i,6]=corrected
        end

        #resetting values
    xChange= Int16(0)
    yChange= Int16(0)
    zChange= Int16(0)
      end#
    end#orthoAxAmove
    offsetIter-=1
    locIterNumb+=1

  end#while
end)#quote
end #singleThreadScan

"""
using single warp we iteratively scan area  where the original point from which dilatation started and led to this result that we are currently analyzing
referenceArrs= (reducedGoldB,reducedSegmB)
                top 6
                bottom 5
                left 2
                right 1
                anterior 3
                posterior 4
"""
macro singleWarpScan(resList, i, referenceArrs , is)

  return  esc(quote
  #this will enable simple way of parallelization as steps will be just executing only those that modulus is compatible with thread id ...

  operationIndex = UInt16(0)

  offsetIter=$resList[$i,6] #storing iteration number
  locIterNumb=0 #will get incemented every time we move centrally in main axis
  while offsetIter>0
    #first we need to  modify x,y or z depending on the direction given in dir variable of result list
    dir = $resList[$i,5]
    xChange-= (dir==2)*offsetIter
    xChange+= (dir==1)*offsetIter

    yChange-= (dir==4)*offsetIter
    yChange+= (dir==3)*offsetIter

    zChange-= (dir==6)*offsetIter
    zChange+= (dir==5)*offsetIter

    #now we have established position over main axis - what is left is to scan by modifying orthogonal axes

    @unroll for orthoAxAmove in -locIterNumb:locIterNumb#so we need to iterate over one of the
      @unroll for orthoAxBmove in [-locIterNumb+ CUDA.abs(orthoAxAmove),locIterNumb-CUDA.abs(orthoAxAmove)]
        #isXMain,isYMain - true if x or in second case y is main axis - we define main axis on the basis of dir variable
        isXMain = (dir==2 || dir==1)
        isYMain = (dir==3 || dir==4)

        #(dir!=2 && dir!=1) so x axis is not main axis
        xChange+= (!isXMain)*orthoAxAmove
        #(dir==2 || dir==1)*orthoAxBmove - so if x is main axis - it means that y is not  and it will be axis A
        yChange+=(!isXMain && !isYMain)*orthoAxBmove+ (isXMain)*orthoAxAmove
        #zChange can never affect axis A  will affect axis B if it is not main axis - so when either x or y is main axis
        zChange+=(!isXMain && !isYMain)*orthoAxBmove

        #here if we will find that there is true in this position we will assume  that those coordinates are coordinates of
        if(referenceArrs[2-$resList[$i,4]][ $resList[$i,1]+xChange,$resList[$i,2]+yChange,$resList[$i,3]+zChange])
            #now on the basis of new coordinates we need to calculate correction for non isometric voxels
            corrected =getCorrectedDistance(yDimSize,zDimSize, $resList[$i,1]+xChange,$resList[$i,2]+yChange,$resList[$i,3]+zChange)
            locSum+=corrected
            $resList[$i,6]=corrected
        end

        #resetting values
    xChange= Int16(0)
    yChange= Int16(0)
    zChange= Int16(0)
      end#
    end#orthoAxAmove
    offsetIter-=1
    locIterNumb+=1

  end#while
end)#quote
end #singleThreadScan

"""
function invoked in case image has non isometric voxels - this ussually is related to the
  fact that voxel is smaller in z dimension than in other dimensions
  IMPORTANT !!! the correction is not perfect as in dilatation step the lack of isometry is not taken into account source voxel may be diffrent than one we would find below, yet in most cases approximation should be good

  x dimension is generally set as a reference dimension and set to be 1 ; y dimension ussually is the same but in theory do not have to be
  z dimension will be given to this function as units where 1 unit is x dimension of voxel in mm
    yDimSize,zDimSize - given in scale where 1 is x dimension of the voxel and presumed to be uniform in whole image
    xDisp,yDisp,zDisp - displacement in x,y,z axis from the covered voxel to source voxel
"""
function getCorrectedDistance(yDimSize,zDimSize, xDisp,yDisp,zDisp)
  return xDisp+ yDimSize*yDisp+ zDimSize*zDisp
end

end#ResultProcess


text:
### Multi-Lane intersection with fixed obstacles

"""
Container type for the environment parameters
"""
@with_kw mutable struct TInterParams
    # geometry
    x_min::Float64 = -20.0
    x_max::Float64 = 20.0
    y_min::Float64 = -20.0
    inter_x::Float64 = 0.
    inter_y::Float64 = 0.
    inter_width::Float64 = 6.0
    inter_height::Float64 = 1.0
    lane_width::Float64 = 3.0
    nlanes_main::Int64 = 2
    nlanes::Int64 = 1

    stop_line::Float64 = 10. # in m along the ego car lane
    # cars
    speed_limit::Float64 = 10.0 # in m/s
    car_rate::Float64 = 0.3 # probability of a car appearing every time steps
end

"""
Generate a T intersection given a geometry specify by TInterParams
"""
function gen_T_roadway(params::TInterParams)
    #=
    ________________________________________________________
        <---- seg 2     main street <----- seg 1
        -----> seg 3                 -----> seg 4
    _________________            ___________________________
                    | seg . seg  |
                    |  5  .  6   |
                    |     .      |
                    |     .      |
                    |     .      |

    =#
    roadway = Roadway()
    # aliases
    x_min = params.x_min
    x_max = params.x_max
    y_min = params.y_min
    inter_x = params.inter_x
    inter_y = params.inter_y
    inter_width = params.inter_width
    inter_height = params.inter_height
    lane_width = params.lane_width
    nlanes_main = params.nlanes_main
    nlanes = params.nlanes

    seg1 = VecSE2(x_max, inter_y+nlanes_main*lane_width - lane_width/2, pi)
    length1 = x_max - inter_x - inter_width
    add_line!(seg1, nlanes_main, length1, roadway)

    seg2 = VecSE2(inter_x - inter_width, inter_y+nlanes_main*lane_width - lane_width/2, pi)
    length2 =  inter_x - inter_width - x_min
    add_line!(seg2, nlanes_main, length2, roadway)

    seg3 = VecSE2(x_min, inter_y-nlanes_main*lane_width + lane_width/2, 0.0)
    length3 = inter_x - x_min - inter_width
    add_line!(seg3, nlanes_main, length3, roadway)

    seg4 = VecSE2(inter_x + inter_width, inter_y-nlanes_main*lane_width + lane_width/2, 0.0)
    length4 = x_max - inter_x - inter_width
    add_line!(seg4, nlanes_main, length4, roadway)

    seg5 = VecSE2(inter_x-nlanes*lane_width + lane_width/2, inter_y - nlanes_main*lane_width - inter_height, -pi/2)
    length5 = inter_y - nlanes_main*lane_width - y_min - inter_height
    add_line!(seg5, nlanes, length5, roadway)

    seg6 = VecSE2(inter_x+nlanes*lane_width - lane_width/2, y_min, pi/2)
    length6 = inter_y - nlanes_main*lane_width - y_min - inter_height
    add_line!(seg6, nlanes, length6, roadway)

    connections = [Connection(1,2), Connection(3,4)]
    push!(connections, Connection(1, 5, 0, [(nlanes_main, 1)]))
    push!(connections, Connection(3, 5, 0, [(1, 1)]))
    # what follows is the ego car left turn
    push!(connections, Connection(6, 2, 0, [(1, nlanes_main)]))
    push!(connections, Connection(6, 4, 0, [(1, 1)]))

    add_junction!(Junction(connections), roadway)
    return roadway
end

function build_obstacles(params::TInterParams)#TODO
    # HEIGHT = 6
    # WIDTH = 6
    #
    # #Left Building
    # left_building = ConvexPolygon(4)
    # top_right = VecSE2(3*params.lane_width/4, -params.lane_width/2 - 1.5, 0.)
    # left_building.pts = [top_right, top_right + VecSE2(-WIDTH, 0., 0.),
    #                 top_right + VecSE2(-WIDTH, -HEIGHT, 0.), top_right + VecSE2(0., -HEIGHT, 0.)]
    # left_building.npts = 4; # a rectangle
    #
    # #Right Building
    # right_building = ConvexPolygon(4)
    # top_right = top_right + VecSE2(top_right.x + 2*params.lane_width + WIDTH, 0., 0.)
    # right_building.pts = [top_right, top_right + VecSE2(-WIDTH, 0., 0.),
    #                 top_right + VecSE2(-WIDTH, -HEIGHT, 0.), top_right + VecSE2(0., -HEIGHT, 0.)]
    # right_building.npts = 4; # a rectangle
    # return [left_building, right_building]
    return Vector{ConvexPolygon}[]
end

"""
T intersection environment
"""
mutable struct IntersectionEnv <: OccludedEnv
    roadway::Roadway
    obstacles::Vector{ConvexPolygon}
    params::TInterParams
    lane_map::Dict{String, Lane}
    priorities::Dict{Tuple{LaneTag, LaneTag}, Bool}
    directions::Dict{Symbol, Tuple{LaneTag, LaneTag}}
end

function IntersectionEnv(params::TInterParams = TInterParams())
    roadway = gen_T_roadway(params)
    lane_map = Dict{String, Lane}()
    obstacles = build_obstacles(params)
    priorities, directions = priority_map(roadway)
    return IntersectionEnv(roadway, obstacles, params, lane_map, priorities, directions)
end

# return a mapping between a label and a route for a 2 lanes T intersection
function priority_map(roadway::Roadway)
    priorities = Dict{Tuple{LaneTag, LaneTag}, Bool}()
    directions = Dict{Symbol, Tuple{LaneTag, LaneTag}}()

    # straight from right
    route = (LaneTag(1,1),LaneTag(2,1))
    priorities[route] = true
    directions[:straight_from_right] = route
    route = (LaneTag(1,2),LaneTag(2,2))
    priorities[route] = true
    directions[:straight_from_right2] = route

    # straight from left
    route = (LaneTag(3,1),LaneTag(4,1))
    priorities[route] = true
    directions[:straight_from_left] = (LaneTag(1,1),LaneTag(2,1))
    route = (LaneTag(3,2),LaneTag(4,2))
    priorities[route] = true
    directions[:straight_from_left2] = (LaneTag(1,2),LaneTag(2,2))

    # right turns
    route = (LaneTag(3,1), LaneTag(5,1))
    priorities[route] = true
    directions[:right_from_left] = route
    # route = (LaneTag(3,2), LaneTag(5,1))
    # priorities[route] = true
    # directions[:right_from_left2] = route
    route = (LaneTag(6,1), LaneTag(4,1))
    priorities[route] = false
    directions[:right_from_bot] = route
    # route = (LaneTag(6,1), LaneTag(4,2))
    # priorities[route] = false
    # directions[:right_from_bot2] = route

    # left turns
    route = (LaneTag(6,1), LaneTag(2,2))
    priorities[route] = false
    route = (LaneTag(6,1), LaneTag(2,1))
    priorities[route] = false
    route = (LaneTag(1,1), LaneTag(5,1))
    priorities[route] = false
    route = (LaneTag(1,2), LaneTag(5,1))
    priorities[route] = false
    return priorities, directions
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
import Base.size, Base.length

export SeisBlock, set_header!, set_traceheader!, set_fileheader!

mutable struct SeisBlock{DT<:Union{IBMFloat32, Float32}}
    fileheader::FileHeader
    traceheaders::AbstractArray{BinaryTraceHeader, 1}
    data::AbstractArray{DT,2}
end

size(block::SeisBlock) = size(block.data)
length(block::SeisBlock) = length(block.traceheaders)

function SeisBlock(data::Array{DT,2}) where {DT<:Union{Float32, IBMFloat32}}

    # Construct FileHeader
    ns, ntraces = size(data)
    fh = FileHeader()
    fh.bfh.ns = ns
    DT==Float32 ? fh.bfh.DataSampleFormat=5 : fh.bfh.DataSampleFormat=1

    return SeisBlock(data, fh)
end

#function SeisBlock(data::Array{DT,2}, fh::FileHeader) where {DT<:Union{Float32, IBMFloat32}}

function SeisBlock(data::Array{DT,2} where {DT<:Union{Float32, IBMFloat32}},
                   fh::FileHeader)

    ns, ntraces = size(data)

    # TODO: Override values in fh.bfh to match data?
#    fh.bfh.ns = ns
#    DT==Float32 ? fh.bfh.DataSampleFormat=5 : fh.bfh.DataSampleFormat=1

    # Construct TraceHeaders
    traceheaders = [BinaryTraceHeader() for i in 1:ntraces]
    set_traceheader!(traceheaders, :ns, ns*ones(Int16, ntraces))

    # Construct Block
    block = SeisBlock(fh, traceheaders, data)

    return block
end


text:
using GenieAuthorisation
using Test

@testset "GenieAuthorisation.jl" begin
    # Write your tests here.
end


text:
# SPDX-License-Identifier: MIT

"""
    abstract type AbstractModel

Abstract type for all (Bayesian) models.
"""
abstract type AbstractModel <: BMObject end

"""
    abstract type AbstractLatentVariableModel

Abstract type for model with latent variable (mixture, pca, ...).
"""
abstract type AbstractLatentVariableModel <: AbstractModel end

"""
    loglikelihood(model::AbstractModel, X)

Return the log-likelihood of `model` of each data sample of `X`.
"""
loglikelihood

"""
    posterior(model::AbstractLatentVariableModel, X)

Return the posterior over the latent variables of `model`.
"""
posterior

"""
    predict(model::AbstractLatentVariableModel, X)

Return the most likely latent variables of `model` given `X`.
"""
predict

"""
    todict(model)

Return a dictionary of the state of the model.
"""
function todict(m::AbstractModel)
    d = Dict()
    names = fieldnames(typeof(m))
    d[:length] = length(names)
    for (i,name) in enumerate(names)
        prop = getproperty(m, name)
        d[Symbol("prop_$(i)_type")] = typeof(prop)
        d[Symbol("prop_$i")] = todict(prop)
    end
    d
end

"""
    fromdict(T, dict)

Return a model of type `T` from the state stored in `dict`.
"""
function fromdict(T::Type{<:AbstractModel}, d::AbstractDict)
    params = []
    for i in 1:d[:length]
        ptype = d[Symbol("prop_$(i)_type")]
        pstate = d[Symbol("prop_$(i)")]
        push!(params, fromdict(ptype, pstate))
    end
    T(params...)
end

"""
    ModelList{N,T}(m1, m2, m3, ...)

Store a list of models for an object's attribute. `N` is the number of
models in the list and `T` is the type of the models. The list is
immutable.
"""
const ModelList{N,T<:AbstractModel} = NTuple{N,T}

max_stars_repo_name: lucasondel/PPCA

max_issues_repo_name: lucasondel/PPCA

max_forks_repo_name: lucasondel/PPCA


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module TrajSaveOnce

export OnceStudy

import Compat.ASCIIString

using AdaptiveStressTesting
using SISLES.GenerativeModel

using CPUTime
using RLESUtils, Obj2Dict, RunCases

import ..DefineSave.trajSave

type OnceStudy
  fileroot::ASCIIString
end

function OnceStudy(;
                   fileroot::AbstractString = "trajSaveONCE")
  OnceStudy(fileroot)
end

function trajSave(study_params::OnceStudy,
                       cases::Cases = Cases(Case());
                       outdir::AbstractString = "./", postproc::Function=identity)

  pmap(case -> begin
         starttime_us = CPUtime_us()
         startnow = string(now())

         sim_params = extract_params!(defineSimParams(), case, "sim_params")
         ast_params = extract_params!(defineASTParams(), case, "ast_params")
         study_params = extract_params!(study_params, case, "study_params")

         sim = defineSim(sim_params)
         ast = defineAST(sim, ast_params)

         simLog = SimLog()
         addObservers!(simLog, ast)

         reward, action_seq = sample(ast)

         notifyObserver(sim, "run_info", Any[reward, sim.md_time, sim.hmd, sim.vmd, sim.label_as_nmac])

         compute_info = ComputeInfo(startnow,
                                    string(now()),
                                    gethostname(),
                                    (CPUtime_us() - starttime_us) / 1e6)

         sav = SaveDict()
         sav["run_type"] = "ONCE"
         sav["compute_info"] = Obj2Dict.to_dict(compute_info)
         sav["sim_params"] = Obj2Dict.to_dict(sim_params)
         sav["ast_params"] = Obj2Dict.to_dict(ast_params)
         sav["sim_log"] = simLog

         fileroot_ = "$(study_params.fileroot)_$(sim.string_id)"
         outfileroot = joinpath(outdir, fileroot_)
         outfile = trajSave(outfileroot, sav)

         #callback for postprocessing
         postproc(outfile)

         return reward
       end,

       cases)
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
#WireInGlobalSettleSelect =

function setGlobalSettlePolicy(fpga,settleWholeHeadstageA,settleWholeHeadstageB, settleWholeHeadstageC,
                                             settleWholeHeadstageD, settleAllHeadstages)

    value = (settleAllHeadstages ? 16 : 0) + (settleWholeHeadstageA ? 1 : 0) + (settleWholeHeadstageB ? 2 : 0) +
            (settleWholeHeadstageC ? 4 : 0) + (settleWholeHeadstageD ? 8 : 0);

    SetWireInValue(fpga,WireInGlobalSettleSelect, value, 0x001f);
    UpdateWireIns(fpga);
}


text:
using Revise
includet(joinpath("..", "init.jl"))

# NOTE: This script was used to make it easier to write the start and end times
# to the temporal blocks. These were written in manually in the database,
# though I suppose it would have been possible to import these parameters

# Load the representative periods and their weights
df = CSV.read(abspath(dataPath, "decision_32_days.csv"))
u = Array(df[:,:periods])[:]
w = Array(df[:,:weights])[:]

# Write the selected day, it's weight plus start and end times to .csv
start = DateTime(2017)
df = DataFrame(
    period = u,
    weight = w,
    start_time = start .+ Day.(u) .- Day(1),
    end_time = start .+ Day.(u)
)

CSV.write(joinpath(dataPath, "32_rep_days_for_Spine.csv"), df)


text:
"""
    po2ai(idx::Vector{Int64}, n::Int64) -> ai_idx::Vector{Int64}

Generator-Notation: Convert PO notation array (padding ones indices) to AI notation array (anti-identity row/column-indices).
"""
function po2ai(idx::Vector{Int64}, n::Int64)::Vector{Int64}
    ai_idx = collect(1:n)
    setdiff!(ai_idx, idx)
    return ai_idx
end

"""
    ai2po(idx::Vector{Int64}, n::Int64) -> po_idx::Vector{Int64}

Generator-Notation: Convert AI notation array (anti-identity row/column-indices) to PO notation array (padding ones indices).
Not really used much but for completenes.

Equivalent to [`po2ai`](@ref).
"""
function ai2po(idx::Vector{Int64}, n::Int64)::Vector{Int64}
    po_idx = collect(1:n)
    setdiff!(po_idx, idx)
    return po_idx
end

"""
    po2matrix(idx::Vector{Int64}, n::Int64) -> SparseMatrixCSC{ℤ₂, Int64}

Generator-Notation: Convert PO to a `n × n` matrix with the anti-identity sub-block at not idx positions
"""
function po2matrix(idx::Vector{Int64}, n::Int64)::SparseMatrixCSC{ℤ₂, Int64}
    k = n - size(idx, 1)
    M = sparse(zeros(ℤ₂, n, n))

    ai_idx = po2ai(idx, n)

    # set the anti-identity subblock
    M[ai_idx, ai_idx] = anti_id(k)
    # set padding ones at diagonal
    for i in idx
        M[i,i] = one(ℤ₂)
    end
    return M
end

"""
    po2latex(idx::Vector{Int64}, n::Int64; po_notation::Bool = false, ai_notation::Bool = false, variable = false, hl_po = "red!10", hl_ai = "green!10", wrapper = "equation*", left = "\\\\left(", right = "\\\\right)") -> String

LaTeX: convert PO to a `n × n` LaTeX-Array with options.

- `idx` (required): Vector with indices in padding ones notation.
- `n` (required): outer matrix size (`n × n`).
- `po_notation = false`: prints out the PO notation if true.
- `ai_notation = false`: print out the AI notation if true.
- `variable = false`: String or false. Prints out the String value if not false.
- `hl_po = "red!10"`: false or LaTeX color string. Highlight padding ones with a given background-color.
- `hl_ai = "green!10`: false or LaTeX color string. Highlight anti-identity with a given background-color.
- `wrapper = "equation"`: LaTeX environment, "equation" or "equation*".
- `left = "\\\\left(`: Left bracket.
- `right = "\\\\right)"`: Right bracket.
"""
function po2latex(
        idx::Vector{Int64}, # PO indices
        n::Int64;
        po_notation::Bool = false, # print PO notation before matrix
        ai_notation::Bool = false, # print AI notation before matrix (and after PO notation)
        variable = false,
        hl_po = "red!10", # (false or LaTeX color string) highlight padding ones with background-color
        hl_ai = "green!10", # (false or LaTeX color string)
        wrapper = "equation*", # LaTeX environment (equation or equation*)
        left = "\\left(",
        right = "\\right)"
    )
    k = n - size(idx, 1)
    M = ones(String, n, n)
    ai_idx = po2ai(idx, n)

    output = "
\\begin{$wrapper}"
    # variable
    if variable != false
        output *= "
    $variable \\equiv "
    end
    # PO notation
    if po_notation == true
        if length(idx) > 1
            output *= "
    $(tuple(idx...))_{PO}^{(n=$n)} = "
        else
            output *= "
    ($(idx...))_{PO}^{(n=$n)} = "
        end
    end
    # AI notation
    if ai_notation == true
        output *= "
    $(tuple(ai_idx...))_{AI}^{(n=$n)} = "
    end

    # set anti-identity subblock
    if hl_ai != false
        hl_ai = "\\cellcolor{$hl_ai}"
    else
        hl_ai = ""
    end
    antiid = fill("$hl_ai 1", k, k)
    for i = 1:k
        antiid[i,i] = "$hl_ai 0"
    end

    M[ai_idx,ai_idx] = antiid

    # set padding ones at diagonal
    if hl_po != false
        hl_po = "\\cellcolor{$hl_po}"
    else
        hl_po = ""
    end
    for i in idx
        M[i,i] = "$hl_po 1"
    end

    # append M to output
    output *= "
    $left\\begin{array}{"
    for i in 1:n
        output *= "c"
    end
    output *= "}"

    for i in 1:size(M,1)
        output *= "
        "
        for j in 1:(size(M,2)-1)
            output *= M[i,j] * " & "
        end
        output *= M[i, end] * " \\\\"
    end

    output *= "
    \\end{array} $right"
    output *= "
\\end{$wrapper}
    "
end

"""
    all_pos2latex(all_pos, n) -> String

LaTeX: Print LaTeX-Arrays for all given base matrices in PO notation in all_pos with a set of default parameters.

See [`po2latex`](@ref) for more information about the possible parameters.
"""
function all_pos2latex(all_pos, n)
    for idx in all_pos
        println(po2latex(idx, n; po_notation=true, hl_po = false))
    end
end

"""
    inv(A::SparseMatrixCSC{ℤ₂})::Matrix{ℤ₂}

Wrapper for the invariant, converts the SparseMatrix back to a matrix, then calls [`Zd_Arithmetics.inv`](https://fhoeddinghaus.github.io/Zd_Arithmetics.jl/operations/#Base.inv-Union{Tuple{Array{%E2%84%A4d{T,%20d},%202}},%20Tuple{d},%20Tuple{T}}%20where%20{T%3C:Integer,%20d}).
"""
inv(A::SparseMatrixCSC{ℤ₂})::Matrix{ℤ₂} = inv(Matrix{ℤ₂}(A))


text:
using JuMP
using Test

const MA = JuMP._MA

include(joinpath(@__DIR__, "utilities.jl"))

@static if !(:JuMPExtension in names(Main))
    include(joinpath(@__DIR__, "JuMPExtension.jl"))
end

# For "expression^3 and unary*"
struct PowVariable <: JuMP.AbstractVariableRef
    pow::Int
end
Base.:^(x::PowVariable, i::Int) = PowVariable(x.pow * i)
Base.:*(x::PowVariable, y::PowVariable) = PowVariable(x.pow + y.pow)
Base.copy(x::PowVariable) = x

function expressions_test(
    ModelType::Type{<:JuMP.AbstractModel},
    VariableRefType::Type{<:JuMP.AbstractVariableRef},
)
    AffExprType = JuMP.GenericAffExpr{Float64,VariableRefType}
    QuadExprType = JuMP.GenericQuadExpr{Float64,VariableRefType}

    @testset "isequal(::GenericAffExpr)" begin
        m = ModelType()
        @variable(m, x)
        @test isequal(x + 1, x + 1)
    end

    @testset "hash(::GenericAffExpr)" begin
        m = ModelType()
        @variable(m, x)
        @test hash(x + 1) == hash(x + 1)
    end

    @testset "drop_zeros!(::GenericAffExpr)" begin
        m = ModelType()
        @variable(m, x[1:2])
        expr = x[1] + x[2] - x[2] + 1
        @test !isequal(expr, x[1] + 1)
        JuMP.drop_zeros!(expr)
        @test isequal(expr, x[1] + 1)
    end

    @testset "iszero(::GenericAffExpr)" begin
        m = ModelType()
        @variable(m, x)
        @test !iszero(x + 1)
        @test !iszero(x + 0)
        @test iszero(0 * x + 0)
        @test iszero(x - x)
    end

    @testset "isequal(::GenericQuadExpr)" begin
        m = ModelType()
        @variable(m, x)
        @test isequal(x^2 + 1, x^2 + 1)
    end

    @testset "hash(::GenericQuadExpr)" begin
        m = ModelType()
        @variable(m, x)
        @test hash(x^2 + 1) == hash(x^2 + 1)
    end

    @testset "drop_zeros!(::GenericQuadExpr)" begin
        m = ModelType()
        @variable(m, x[1:2])
        expr = x[1]^2 + x[2]^2 - x[2]^2 + x[1] + x[2] - x[2] + 1
        @test !isequal(expr, x[1]^2 + x[1] + 1)
        JuMP.drop_zeros!(expr)
        @test isequal(expr, x[1]^2 + x[1] + 1)
    end

    @testset "iszero(::GenericQuadExpr)" begin
        m = ModelType()
        @variable(m, x)
        @test !iszero(x^2 + 1)
        @test !iszero(x^2 + 0)
        @test !iszero(x^2 + 0 * x + 0)
        @test iszero(0 * x^2 + 0 * x + 0)
        @test iszero(x^2 - x^2)
    end

    @testset "value for GenericAffExpr" begin
        expr1 = JuMP.GenericAffExpr(3.0, 3 => -5.0, 2 => 4.0)
        @test @inferred(JuMP.value(expr1, -)) == 10.0
        expr2 = JuMP.GenericAffExpr{Int,Int}(2)
        @test typeof(@inferred(JuMP.value(expr2, i -> 1.0))) == Float64
        @test @inferred(JuMP.value(expr2, i -> 1.0)) == 2.0
    end

    @testset "value for GenericQuadExpr" begin
        # 1 + 2x(1) + 3x(2)
        affine_term = JuMP.GenericAffExpr(1.0, 1 => 2.0, 2 => 3.0)
        # 1 + 2x(1) + 3x(2) + 4x(1)^2 + 5x(1)*x(2) + 6x(2)^2
        expr = JuMP.GenericQuadExpr(
            affine_term,
            JuMP.UnorderedPair(1, 1) => 4.0,
            JuMP.UnorderedPair(1, 2) => 5.0,
            JuMP.UnorderedPair(2, 2) => 6.0,
        )
        @test typeof(@inferred(JuMP.value(expr, i -> 1.0))) == Float64
        @test @inferred(JuMP.value(expr, i -> 1.0)) == 21
        @test @inferred(JuMP.value(expr, i -> 2.0)) == 71
    end

    @testset "add_to_expression!(::GenericAffExpr{C,V}, ::V)" begin
        aff = JuMP.GenericAffExpr(1.0, :a => 2.0)
        @test JuMP.isequal_canonical(
            JuMP.add_to_expression!(aff, :b),
            JuMP.GenericAffExpr(1.0, :a => 2.0, :b => 1.0),
        )
    end

    @testset "add_to_expression!(::GenericAffExpr{C,V}, ::C)" begin
        aff = JuMP.GenericAffExpr(1.0, :a => 2.0)
        @test JuMP.isequal_canonical(
            JuMP.add_to_expression!(aff, 1.0),
            JuMP.GenericAffExpr(2.0, :a => 2.0),
        )
    end

    @testset "linear_terms(::AffExpr)" begin
        m = ModelType()
        @variable(m, x[1:10])

        aff = 1 * x[1] + 2 * x[2]
        k = 0
        @test length(linear_terms(aff)) == 2
        for (coeff, var) in linear_terms(aff)
            if k == 0
                @test coeff == 1
                @test var === x[1]
            elseif k == 1
                @test coeff == 2
                @test var === x[2]
            end
            k += 1
        end
        @test k == 2
    end

    @testset "linear_terms(::AffExpr) for empty expression" begin
        k = 0
        aff = zero(AffExprType)
        @test length(linear_terms(aff)) == 0
        for (coeff, var) in linear_terms(aff)
            k += 1
        end
        @test k == 0
    end

    @testset "coefficient(aff::AffExpr, v::VariableRef)" begin
        m = ModelType()
        x = @variable(m, x)
        y = @variable(m, y)
        aff = @expression(m, 1.0 * x)
        @test coefficient(aff, x) == 1.0
        @test coefficient(aff, y) == 0.0
    end

    @testset "coefficient(aff::AffExpr, v1::VariableRef, v2::VariableRef)" begin
        m = ModelType()
        x = @variable(m, x)
        aff = @expression(m, 1.0 * x)
        @test coefficient(aff, x, x) == 0.0
    end

    @testset "coefficient(quad::QuadExpr, v::VariableRef)" begin
        m = ModelType()
        x = @variable(m, x)
        y = @variable(m, y)
        z = @variable(m, z)
        quad = @expression(m, 6.0 * x^2 + 5.0 * x * y + 2.0 * y + 3.0 * x)
        @test coefficient(quad, x) == 3.0
        @test coefficient(quad, y) == 2.0
        @test coefficient(quad, z) == 0.0
    end

    @testset "coefficient(quad::Quad, v1::VariableRef, v2::VariableRef)" begin
        m = ModelType()
        x = @variable(m, x)
        y = @variable(m, y)
        z = @variable(m, z)
        quad = @expression(m, 6.0 * x^2 + 5.0 * x * y + 2.0 * y + 3.0 * x)
        @test coefficient(quad, x, y) == 5.0
        @test coefficient(quad, x, x) == 6.0
        @test coefficient(quad, x, y) == coefficient(quad, y, x)
        @test coefficient(quad, z, z) == 0.0
    end

    @testset "MA.add_mul!(ex::Number, c::Number, x::GenericAffExpr)" begin
        aff = MA.add_mul!(1.0, 2.0, JuMP.GenericAffExpr(1.0, :a => 1.0))
        @test JuMP.isequal_canonical(aff, JuMP.GenericAffExpr(3.0, :a => 2.0))
    end

    @testset "MA.add_mul!(ex::Number, c::Number, x::GenericQuadExpr) with c == 0" begin
        quad = MA.add_mul!(2.0, 0.0, QuadExprType())
        @test JuMP.isequal_canonical(quad, convert(QuadExprType, 2.0))
    end

    @testset "MA.add_mul!(ex::Number, c::VariableRef, x::VariableRef)" begin
        model = ModelType()
        @variable(model, x)
        @variable(model, y)
        @test_expression_with_string MA.add_mul(5.0, x, y) "x*y + 5"
        @test_expression_with_string MA.add_mul!(5.0, x, y) "x*y + 5"
    end

    @testset "MA.add_mul!(ex::Number, c::T, x::T) where T<:GenericAffExpr" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(1.0, 2x, x + 1) "2 x² + 2 x + 1"
        @test_expression_with_string MA.add_mul!(1.0, 2x, x + 1) "2 x² + 2 x + 1"
    end

    @testset "MA.add_mul!(ex::Number, c::GenericAffExpr{C,V}, x::V) where {C,V}" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(1.0, 2x, x) "2 x² + 1"
        @test_expression_with_string MA.add_mul!(1.0, 2x, x) "2 x² + 1"
    end

    @testset "MA.add_mul!(ex::Number, c::GenericQuadExpr, x::Number)" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(0.0, x^2, 1.0) "x²"
        @test_expression_with_string MA.add_mul!(0.0, x^2, 1.0) "x²"
    end

    @testset "MA.add_mul!(ex::Number, c::GenericQuadExpr, x::Number) with c == 0" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(0.0, x^2, 0.0) "0"
        @test_expression_with_string MA.add_mul!(0.0, x^2, 0.0) "0"
    end

    @testset "MA.add_mul!(aff::AffExpr,c::VariableRef,x::AffExpr)" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(2x, x, x + 1) "x² + 3 x"
        @test_expression_with_string MA.add_mul!(2x, x, x + 1) "x² + 3 x"
    end

    @testset "MA.add_mul!(aff::GenericAffExpr{C,V},c::GenericAffExpr{C,V},x::Number) where {C,V}" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(2x, x, 1) "3 x"
        @test_expression_with_string MA.add_mul!(2x, x, 1) "3 x"
    end

    @testset "MA.add_mul!(aff::GenericAffExpr{C,V}, c::GenericQuadExpr{C,V}, x::Number) where {C,V}" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(2x, x^2, 1) "x² + 2 x"
        @test_expression_with_string MA.add_mul!(2x, x^2, 1) "x² + 2 x"
    end

    @testset "MA.add_mul!(aff::GenericAffExpr{C,V}, c::GenericQuadExpr{C,V}, x::Number) where {C,V} with x == 0" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(2x, x^2, 0) "2 x"
        @test_expression_with_string MA.add_mul!(2x, x^2, 0) "2 x"
    end

    @testset "MA.add_mul!(aff::GenericAffExpr{C,V}, c::Number, x::GenericQuadExpr{C,V}) where {C,V} with c == 0" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(2x, 0, x^2) "2 x"
        @test_expression_with_string MA.add_mul!(2x, 0, x^2) "2 x"
    end

    @testset "MA.add_mul!(ex::GenericAffExpr{C,V}, c::GenericAffExpr{C,V}, x::GenericAffExpr{C,V}) where {C,V}" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(2x, x + 1, x + 0) "x² + 3 x"
        @test_expression_with_string MA.add_mul!(2x, x + 1, x + 0) "x² + 3 x"
    end

    @testset "MA.add_mul!(quad::GenericQuadExpr{C,V},c::GenericAffExpr{C,V},x::Number) where {C,V}" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(x^2, x + 1, 1) "x² + x + 1"
        @test_expression_with_string MA.add_mul!(x^2, x + 1, 1) "x² + x + 1"
    end

    @testset "MA.add_mul!(quad::GenericQuadExpr{C,V},c::V,x::GenericAffExpr{C,V}) where {C,V}" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(x^2, x, x + 1) "2 x² + x"
        @test_expression_with_string MA.add_mul!(x^2, x, x + 1) "2 x² + x"
    end

    @testset "MA.add_mul!(quad::GenericQuadExpr{C,V},c::GenericQuadExpr{C,V},x::Number) where {C,V}" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(x^2 + x, x^2 + x, 2.0) "3 x² + 3 x"
        @test_expression_with_string MA.add_mul!(x^2 + x, x^2 + x, 2.0) "3 x² + 3 x"
    end

    @testset "MA.add_mul!(ex::GenericQuadExpr{C,V}, c::GenericAffExpr{C,V}, x::GenericAffExpr{C,V}) where {C,V}" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string MA.add_mul(x^2 + x, x + 0, x + 1) "2 x² + 2 x"
        @test_expression_with_string MA.add_mul!(x^2 + x, x + 0, x + 1) "2 x² + 2 x"
    end

    @testset "(+)(::AffExpr)" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string (+)(x + 1) "x + 1"
    end

    @testset "(+)(::QuadExpr)" begin
        model = ModelType()
        @variable(model, x)
        @test_expression_with_string (+)(x^2 + 1) "x² + 1"
    end

    @testset "sum(::Vector{VariableRef})" begin
        model = ModelType()
        @variable(model, x[1:2])
        @test_expression_with_string sum(x) "x[1] + x[2]"
    end

    @testset "expression^3 and unary*" begin
        model = ModelType()
        x = PowVariable(1)
        # Calls (*)((x*x)^6)
        y = @expression model (x * x)^3
        @test y.pow == 6
        z = @inferred (x * x)^3
        @test z.pow == 6
    end

    @testset "ndims(::QuadExpr)" begin
        model = ModelType()
        @variable(model, x)
        @test ndims(x^2 + 1) == 0
    end
end

@testset "Expressions for JuMP.Model" begin
    expressions_test(Model, VariableRef)
end

@testset "Expressions for JuMPExtension.MyModel" begin
    expressions_test(JuMPExtension.MyModel, JuMPExtension.MyVariableRef)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using BinDeps
@static if VERSION <= v"0.7-"
	import Compat: Sys
end

@BinDeps.setup

const libui = library_dependency("libui")
const libuiVer = "4"
const libuiFilebase = "alpha$(libuiVer)"

const prefix = joinpath(BinDeps.depsdir(libui), "usr")

provides(Sources,
         URI("http://w3.ualg.pt/%7Ejluis/t/$(libuiFilebase).tar.gz"),
         libui)

const prefix = joinpath(BinDeps.depsdir(libui), "usr")

const srcpath = if Sys.islinux()
    joinpath(BinDeps.srcdir(libui), libuiFilebase, "linux_amd64", "shared", "libui.so.0")
elseif Sys.isapple()
    joinpath(BinDeps.srcdir(libui), libuiFilebase, "darwin_amd64", "shared", "libui.A.dylib")
elseif Sys.iswindows() #TODO: Check this on a CI service
    joinpath(BinDeps.srcdir(libui), libuiFilebase, "windows_amd64", "shared", "libui.dll")
else
    error("Platform not supported.")
end

const binaryNameTarget = if Sys.islinux()
	"libui.so"
elseif Sys.isapple()
	"libui.dylib"
elseif Sys.iswindows()
	"libui.dll"
else
	error("Platform not supported.")
end

provides(SimpleBuild,
		 (@build_steps begin
		  GetSources(libui)
		  CreateDirectory(joinpath(prefix, "lib"))
		  `cp -v $(srcpath) $(joinpath(prefix, "lib", binaryNameTarget))`
		  end
		  ), libui)

@BinDeps.install Dict(:libui => :libui)


text:
@testset "trace" begin
    @testset "basic usage" begin
        global_tracer_provider!(
            TracerProvider(
                span_processor = CompositSpanProcessor(
                    SimpleSpanProcessor(InMemoryExporter()),
                ),
            ),
        )

        tracer = Tracer()

        with_span("test_async", tracer) do
            @sync for i in 1:5
                @async with_span("asyn_span_$i", tracer) do
                    current_span()["my_id"] = i
                end
            end
        end

        p = global_tracer_provider()
        sp = p.span_processor.span_processors[1]
        @test length(sp.span_exporter.pool) == 6

        push!(p, SimpleSpanProcessor(ConsoleExporter()))

        @test_throws ErrorException with_span("test", tracer) do
            span_name!("TEST")  # one can change the name of current span
            println("hello world from $(span_name()) [$(span_status())]!")
            try
                throw(ErrorException("???"))
            catch e
                push!(current_span(), e; is_rethrow_followed = true)
                rethrow(e)
            end
        end

        force_flush!(p)
        shut_down!(p)
        @test length(sp.span_exporter.pool) == 0

        with_span("test", tracer) do
            @test span_context() === INVALID_SPAN_CONTEXT
        end
    end

    @testset "sampling" begin
        p = TracerProvider(
            span_processor = SimpleSpanProcessor(InMemoryExporter()),
            sampler = ALWAYS_OFF,
        )

        t = Tracer(; provider = p)

        with_span("test", t) do
            @test span_context() === INVALID_SPAN_CONTEXT
        end

        p = TracerProvider(
            span_processor = SimpleSpanProcessor(InMemoryExporter()),
            sampler = DEFAULT_OFF,
        )

        t = Tracer(; provider = p)

        with_span("test", t) do
            @test span_context() === INVALID_SPAN_CONTEXT
        end

        p = TracerProvider(
            span_processor = SimpleSpanProcessor(InMemoryExporter()),
            sampler = TraceIdRatioBased(0.5),
            id_generator = RandomIdGenerator(; rng = MersenneTwister(123)),
        )

        t = Tracer(; provider = p)
        n_no_op_spans = 0
        for _ in 1:1_000
            with_span("test", t) do
                if span_context() === INVALID_SPAN_CONTEXT
                    n_no_op_spans += 1
                end
            end
        end
        @test isapprox(n_no_op_spans / 1_000, 0.5; atol = 0.1)
    end
end


text:
# This file contains types and structures for productivity DEA model
"""
    AbstractProductivityDEAlModel
An abstract type representing a productivity DEA model.
"""
abstract type AbstractProductivityDEAModel <: AbstractDEAModel end

"""
    nperiods(model::AbstractProductivityDEAModel)
Return number of time periods of a productivity DEA model.
# Examples
```jldoctest
julia> X = Array{Float64,3}(undef, 5, 1, 2);

julia> X[:, :, 1] = [2; 3; 5; 4; 4];

julia> X[:, :, 2] = [1; 2; 4; 3; 4];

julia> Y = Array{Float64,3}(undef, 5, 1, 2);

julia> Y[:, :, 1] = [1; 4; 6; 3; 5];

julia> Y[:, :, 2] = [1; 4; 6; 3; 3];

julia> mprod = malmquist(X, Y);

julia> nperiods(mprod)
2
```
"""
nperiods(model::AbstractProductivityDEAModel) = model.periods

"""
    prodchange(model::AbstractProductivityDEAModel)
Return productivity change of a productivity change DEA model.

# Optional Arguments
- `type=Prod`: component of productivity change to return.

Type specification:
- `:Prod`: productivity change.
- `:EC`: efficiency change.
- `:TC`: technological change.

# Examples
```jldoctest
julia> X = Array{Float64,3}(undef, 5, 1, 2);

julia> X[:, :, 1] = [2; 3; 5; 4; 4];

julia> X[:, :, 2] = [1; 2; 4; 3; 4];

julia> Y = Array{Float64,3}(undef, 5, 1, 2);

julia> Y[:, :, 1] = [1; 4; 6; 3; 5];

julia> Y[:, :, 2] = [1; 4; 6; 3; 3];

julia> mprod = malmquist(X, Y);

julia> prodchange(mprod)
×1 Matrix{Float64}:
 2.0
 1.5
 1.25
 1.3333333333333333
 0.6000000000000001

julia> prodchange(mprod, :EC)
5×1 Matrix{Float64}:
 1.3333333333333333
 1.0
 0.8333333333333334
 0.8888888888888888
 0.4

julia> prodchange(mprod, :TC)
5×1 Matrix{Float64}:
 1.5
 1.5
 1.5
 1.5
 1.5
```
"""
function prodchange(model::AbstractProductivityDEAModel, type::Symbol = :Prod)::Matrix

    if type == :Prod
        return model.Prod
    end

    if type == :EC
        if isdefined(model, :EC)
            return model.EC
        end
    end

    if type == :TC
        if isdefined(model, :TC)
            return model.TC
        end
    end

    throw(ArgumentError("$(typeof(model)) has no productivity change component $(type)"));

end


text:
rand_inbounds(r::MRG32k3a, ::Random.CloseOpen12_64) = rand(r)+1.0
rand_inbounds(r::MRG32k3a, ::Random.CloseOpen01_64=Random.CloseOpen01()) = rand(r)

############################
rand_inbounds(r::MRG32k3a, ::Random.UInt52Raw{T}) where {T<:Random.BitInteger} = reinterpret(UInt64, rand_inbounds(r, Random.CloseOpen12())) % T

function rand(r::MRG32k3a, x::Random.SamplerTrivial{Random.UInt52Raw{UInt64}})
    rand_inbounds(r, x[])
end

function rand(r::MRG32k3a, ::Random.SamplerTrivial{Random.UInt2x52Raw{UInt128}})
    rand_inbounds(r, Random.UInt52Raw(UInt128)) << 64 | rand_inbounds(r, Random.UInt52Raw(UInt128))
end

function rand(r::MRG32k3a, ::Random.SamplerTrivial{Random.UInt104Raw{UInt128}})
    reserve(r, 2)
	xor(rand_inbounds(r, Random.UInt52Raw(UInt128)) << 52, rand_inbounds(r, Random.UInt52Raw(UInt128)))
end

#### floats

rand(r::MRG32k3a, sp::Random.SamplerTrivial{Random.CloseOpen12_64}) = rand_inbounds(r, sp[])


text:
using NeRCA
using Test

const DETX = joinpath(@__DIR__, "data", "detx_v3.detx")
const H5FILE = joinpath(@__DIR__, "data", "mupage.h5")

@testset "calibration" begin
    calib = NeRCA.read_calibration(DETX)

    @test 23 == calib.det_id
    @test 6 == length(values(calib.pos))
    @test 6 == length(values(calib.dir))
    @test 6 == length(values(calib.t0))
    @test 6 == length(values(calib.du))
    @test 6 == length(values(calib.floor))
    @test 18 == length(values(calib.omkeys))
    @test 23.9 ≈ calib.max_z
    @test 2 == calib.n_dus
    @test 3.6 ≈ calib.pos[3][2].z
end

@testset "reading KM3HDF5 hits" begin
    hits = read_hits(H5FILE, 0)
    @test 4886 == length(hits)
    @test 3.3263139e7 ≈ hits[end].t
    hits = read_hits(H5FILE, 2)
    @test 4219 == length(hits)
    @test 5.1047085e7 ≈ hits[end].t
end

@testset "reading KM3HDF5 MC hits" begin
    hits = NeRCA.read_mchits(H5FILE, 0)
    @test 399 == length(hits)
    @test 4429.716844306684 ≈ hits[end].t
end

@testset "reading KM3HDF5 MC tracks" begin
    mc_tracks = NeRCA.read_mctracks(H5FILE)
    @test 3 == length(mc_tracks)
    @test 11 == length(mc_tracks[0])
    @test 3 == length(mc_tracks[2])
end


text:
using IncrementalInference
using IncrementalInference.Optim
using Test

##
@testset "test parametric mixture" begin

## Test simple mixture prior

fg = LocalDFG(;solverParams=SolverParams(algorithms=[:default, :parametric]))

addVariable!(fg, :x0, ContinuousScalar)

p = Mixture(Prior(I), [Normal(0.8, 0.4), Normal(1.0, 0.1)], Categorical([0.1; 0.9]))
f = addFactor!(fg, [:x0], p)

addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0,:x1], LinearRelative(Normal(1.0,0.1)))

options = Optim.Options(time_limit = 100,
                    iterations = 10000,
                    show_trace = true,
                    show_every = 10,
                    allow_f_increases=true,
                    g_tol = 1e-6,
                    )

algorithm = Optim.NelderMead
vardict, result, varIds, Σ = solveGraphParametric(fg; options, algorithm)

@test isapprox(vardict[:x0].val[1], 1, atol = 0.05)
@test isapprox(vardict[:x0].cov[1], 0.01, atol = 0.001)

##
fg = LocalDFG(;solverParams=SolverParams(algorithms=[:default, :parametric]))

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0], Prior(Normal(0.0,0.1)))

mlr = Mixture(LinearRelative(I), [Normal(-1.0, 0.2), Normal(1.0, 0.1)], Categorical([0.5; 0.5]))

addFactor!(fg, [:x0,:x1], mlr)

addVariable!(fg, :l1, ContinuousScalar)
p = Mixture(Prior(I), [Normal(-1.5, 0.1), Normal(0.9, 0.2)], Categorical([0.5; 0.5]))
# p = Prior(Normal(0.9, 0.1))
addFactor!(fg, [:l1], p)

addFactor!(fg, [:x1,:l1], LinearRelative(Normal(0.0,0.1)))

vardict, result, varIds, Σ = solveGraphParametric(fg)
@test isapprox(vardict[:x0].val[1], 0, atol = 0.1)
@test isapprox(vardict[:x1].val[1], 1, atol = 0.1)
@test isapprox(vardict[:l1].val[1], 1, atol = 0.1)

## ContinuousEuclid(2) prior

fg = LocalDFG(;solverParams=SolverParams(algorithms=[:default, :parametric]))

addVariable!(fg, :x0, ContinuousEuclid(2))

p = Mixture(Prior(MvNormal(2,1.0)), [MvNormal([0.8, 0.5], [0.4, 0.4]), MvNormal([1.0, 0.5], [0.1, 0.1])], Categorical([0.1; 0.9]))
f = addFactor!(fg, [:x0], p)

vardict, result, varIds, Σ = solveGraphParametric(fg)
vardict
@test isapprox(vardict[:x0].val[1], 1.0, atol = 0.01)
@test isapprox(vardict[:x0].val[2], 0.5, atol = 0.01)
@test isapprox(vardict[:x0].cov[1,1], 0.01, atol = 0.001)
@test isapprox(vardict[:x0].cov[2,2], 0.01, atol = 0.001)

# x = collect(0.5:0.01:1.5)
# y = collect(0.0:0.01:1.0)
# cfm = IIF.CalcFactorMahalanobis(f)
# xy = Vector{Float64}[]
# x = 0.5:0.01:1.5
# y = 0:0.01:1
# for x=x,y=y
#     push!(xy, [x,y])
# end
# r = cfm.(xy)
# surface(x,y,r)

##
end

if false

fg = LocalDFG(;solverParams=SolverParams(algorithms=[:default, :parametric]))

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0], Prior(Normal(0.0,0.1)))

mlr = Mixture(LinearRelative(I), [Normal(-1.0, 0.2), Normal(1.0, 0.1)], Categorical([0.5; 0.5]))
# mlr = Mixture(LinearRelative(I), [Normal(-0.2, 0.05), Normal(0.2, 0.1)], Categorical([0.5; 0.5]))
# mlr = Mixture(LinearRelative(I), [Normal(0.0, 0.9), Normal(0.0, 0.1)], Categorical([0.5; 0.5]))

# mlr = Mixture(LinearRelative(I), [Normal(-0.2, 0.1), Normal(0.2, 0.1)], Categorical([0.9; 0.1]))
# mlr = Mixture(LinearRelative(I), [Normal(-0.2, 0.15), Normal(0.2, 0.1)], Categorical([0.5; 0.5]))
# mlr = Mixture(LinearRelative(I), [Normal(-1.0, 0.2), Normal(1.0, 0.1)], Categorical([0.5; 0.5]))

addFactor!(fg, [:x0,:x1], mlr)

f = fg[:x0x1f1]
cfm = IIF.CalcFactorMahalanobis(f)
x1 = collect(-2:0.01:2)
# x1 = reverse(collect(-2:0.01:2))
# x1 = collect(-0.5:0.01:0.5)
r1 = cfm.(Ref([0.0]), x1)
# plot!(x1,r)
plot(x1,r1)

##
addVariable!(fg, :l1, ContinuousScalar)
p = Mixture(Prior(I), [Normal(-1.5, 0.1), Normal(0.9, 0.2)], Categorical([0.5; 0.5]))
# p = Prior(Normal(0.9, 0.1))
addFactor!(fg, [:l1], p)

addFactor!(fg, [:x1,:l1], LinearRelative(Normal(0.0,0.1)))

f = fg[:l1f1]
cfm = IIF.CalcFactorMahalanobis(f)
x1 = collect(-2:0.01:2)
r2 = cfm.(x1)
plot!(x1,r2)
# plot!(x1,r1+r2)

##
initAll!(fg)
solveTree!(fg)
plotKDE(fg, :x1)
plotKDE(fg, :l1)
##
using IncrementalInference.Optim

options = Optim.Options(time_limit = 100,
                    iterations = 1000,
                    show_trace = true,
                    show_every = 1,
                    allow_f_increases=true,
                    g_tol = 1e-6,
                    )

algorithm = Optim.NelderMead
algorithm = Optim.BFGS
algorithmkwargs=(linesearch = Optim.LineSearches.HagerZhang(),)
# algorithmkwargs=(linesearch = Optim.LineSearches.Static(),)
vardict, result, varIds, Σ = solveGraphParametric(fg; algorithm, options, algorithmkwargs)
vardict

##

IIF.updateParametricSolution!(fg, vardict)

##

end

##
if false
using RoME

fg = LocalDFG(;solverParams=SolverParams(algorithms=[:default, :parametric]))

pr_noise = [0.01, 0.01, 0.001]
od_noise_1 = [0.5; 0.5; 0.2]
od_noise_2 = [0.05; 0.05; 0.02]

#x0 prior
addVariable!(fg, :x0, Pose2)
prpo = MvNormal([0,0,-pi], pr_noise)
addFactor!(fg, [:x0], PriorPose2(prpo))
# addFactor!(fg, [:x0], PriorPose2(MvNormal(rand(prpo), pr_noise)))

#x0 to x1
addVariable!(fg, :x1, Pose2)
pp1 = MvNormal([1.0,0,0], od_noise_1)
pp2 = MvNormal([1.0,0,0], od_noise_2)
mpp = Mixture(Pose2Pose2(I), [pp1, pp2], Categorical([0.5; 0.5]))

f = addFactor!(fg, [:x0,:x1], mpp)

cfm = IIF.CalcFactorMahalanobis(f)

##
xyθ = Vector{Float64}[]
for x = 0.5:0.01:1.5, y = 0, θ = 0
    push!(xyθ, [x,y,θ])
end

x = 0.5:0.01:1.5

r = cfm.(Ref([0.0; 0.0; 0.0]), xyθ)
plot(x,r)

##
#x1 to x2
# addVariable!(fg, :x2, Pose2)
# pp1 = MvNormal([1.0,0,0], od_noise_1)
# pp2 = MvNormal([1.0,0,0], od_noise_2)
# mpp = Mixture(Pose2Pose2(I), [pp1, pp2], Categorical([0.5; 0.5]))
# addFactor!(fg, [:x1,:x2], mpp)

# vardict, result, varIds, Σ = solveFactorGraphParametric(fg)
# vardict

##

using IncrementalInference.Optim

options = Optim.Options(time_limit = 100,
                    iterations = 10000,
                    show_trace = true,
                    show_every = 1,
                    allow_f_increases=true,
                    g_tol = 1e-6,
                    )

algorithm = Optim.NelderMead
vardict, result, varIds, Σ = solveGraphParametric(fg; algorithm, autodiff=:finite)
vardict

end


text:
module TapirDevUtils

export @ci_tapir,
    @ircode_tapir,
    child_task_info,
    ci_child_tasks,
    ci_tapir,
    ircode_child_tasks,
    ircode_tapir,
    newworld_compiler

using InteractiveUtils:
    InteractiveUtils, code_typed, gen_call_with_extracted_types_and_kwargs

try
    include("unstable.jl")
catch err
    @error "Cannot import unstable API" exception = (err, catch_backtrace())
end

function ci_tapir(f, types)
    @nospecialize
    (ci,), = code_typed(f, types)
    mch = Base._which(Base.tuple_type_cons(typeof(f), types))
    mi = Core.Compiler.specialize_method(mch.method, mch.spec_types, mch.sparams)
    return Core.Compiler.lower_tapir(mi, ci)
end

"""
    @ci_tapir f(args...)

Get `CodeInfo` after processed for lowering to Tapir. This is the actual
`CodeInfo` compiled to LLVM (while `@code_typed` still includes Tapir
constructs.)
"""
macro ci_tapir(args...)
    gen_call_with_extracted_types_and_kwargs(__module__, ci_tapir, args)
end

function ircode_tapir(f, types)
    @nospecialize
    (ci,), = code_typed(f, types)
    mch = Base._which(Base.tuple_type_cons(typeof(f), types))
    mi = Core.Compiler.specialize_method(mch.method, mch.spec_types, mch.sparams)
    return Core.Compiler.lower_tapir_to_ircode(mi, ci)
end

"""
    @ircode_tapir f(args...)

Get `IRCode` after processed for lowering to Tapir.
"""
macro ircode_tapir(args...)
    gen_call_with_extracted_types_and_kwargs(__module__, ircode_tapir, args)
end

"""
    newworld_compiler()

Use Revised `Core.Compiler` for real compilation (execution).

It just calls:

    ccall(:jl_set_typeinf_func, Cvoid, (Any,), Core.Compiler.typeinf_ext_toplevel)
"""
newworld_compiler() =
    ccall(:jl_set_typeinf_func, Cvoid, (Any,), Core.Compiler.typeinf_ext_toplevel)

"""
    UnsafeTapir

A namespace that provides Tapir-like API but with simpler but more "unsafe" code
(e.g., no exception handling).  Useful for debugging IR.

Usage:
    using TapirDevUtils: UnsafeTapir as Tapir
"""
module UnsafeTapir
using Base.Experimental: Tapir

try
    using Base.Experimental.Tapir: @output
catch
end

const tokenname = gensym(:token)

macro spawn(block)
    var = esc(tokenname)
    block = Expr(:block, __source__, block)
    :(Tapir.@spawnin $var $(esc(block)))
end

macro sync(block)
    var = esc(tokenname)
    block = Expr(:block, __source__, block)
    quote
        let $var = Tapir.@syncregion(), ans = $(esc(block))
            Tapir.@sync_end($var)
            ans
        end
    end
end

end # module UnsafeTapir

const Requires = try
    Base.require(Base.PkgId(Base.UUID(0xae029012a4dd51049daad747884805df), "Requires"))
catch
    nothing
end

@static if Requires isa Module
    using .Requires: @require
end

@static if Requires isa Module
    function __init__()
        @require Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f" include("cthulhu.jl")
    end
end

end


text:
export
    PEPSNetwork,
    generate_boundary,
    node_from_index,
    conditional_probability

const IntOrRational = Union{Int, Rational{Int}}

struct PEPSNetwork{F} <: AbstractGibbsNetwork{NTuple{2, Int}, NTuple{2, Int}} where F
    factor_graph::LabelledGraph{T, NTuple{2, Int}} where T
    network_graph::LabelledGraph{S, NTuple{2, Int}} where S
    vertex_map::Function
    m::Int
    n::Int
    nrows::Int
    ncols::Int
    β::Real
    bond_dim::Int
    var_tol::Real
    sweeps::Int
    gauges::Dict{Tuple{Rational, IntOrRational}, Vector{Float64}}
    tensor_spiecies::Dict{NTuple{2, IntOrRational}, Symbol}
    columns_MPO::NTuple{N, IntOrRational} where N
    layers_MPS::NTuple{M, IntOrRational} where M
    layers_left_env::NTuple{K, IntOrRational} where K
    layers_right_env::NTuple{L, IntOrRational} where L

    function PEPSNetwork{F}(
        m::Int,
        n::Int,
        factor_graph::LabelledGraph,
        transformation::LatticeTransformation;
        β::Real,
        bond_dim::Int=typemax(Int),
        var_tol::Real=1E-8,
        sweeps::Int=4,
        columns_MPO = (-1//2, 0),  # from left to right
        layers_MPS=(4//6, 3//6, 0, -1//6),  # from bottom to top
        layers_left_env=(4//6, 3//6),
        layers_right_env=(0, -3//6)
    )  where F
        vmap = vertex_map(transformation, m, n)
        ng = F ? cross_lattice(m, n) : peps_lattice(m, n)
        nrows, ncols = transformation.flips_dimensions ? (n, m) : (m, n)

        if !is_compatible(factor_graph, ng)
            throw(ArgumentError("Factor graph not compatible with given network."))
        end

        _types = (:site, :central_h, :central_v, :gauge_h)
        _gauges = Dict{Tuple{Rational, Int}, Vector{Float64}}()
        _tensor_spiecies = Dict{NTuple{2, IntOrRational}, Symbol}()

        network = new(
            factor_graph,
            ng,
            vmap,
            m,
            n,
            nrows,
            ncols,
            β,
            bond_dim,
            var_tol,
            sweeps,
            _gauges,
            _tensor_spiecies,
            columns_MPO,
            layers_MPS,
            layers_left_env,
            layers_right_env
        )

        assign_tensors!(network)
        update_gauges!(network, :id)
        network
    end
end

function peps_lattice(m::Int, n::Int)
    labels = [(i, j) for j ∈ 1:n for i ∈ 1:m]
    LabelledGraph(labels, grid((m, n)))
end

function cross_lattice(m::Int, n::Int)
    labels = [(i, j) for j ∈ 1:n for i ∈ 1:m]
    lg = LabelledGraph(labels, grid((m, n)))
    for i ∈ 1:m-1, j ∈ 1:n-1
        add_edge!(lg, (i, j), (i+1, j+1))
        add_edge!(lg, (i+1, j), (i, j+1))
    end
    lg
end

function projectors(network::PEPSNetwork{false}, vertex::NTuple{2, Int})
    i, j = vertex
    neighbours = ((i, j-1), (i-1, j), (i, j+1), (i+1, j))
    projector.(Ref(network), Ref(vertex), neighbours)
end

function projectors(network::PEPSNetwork{true}, vertex::NTuple{2, Int})
    i, j = vertex
    neighbours = (
                    (
                        (i+1, j-1), (i, j-1), (i-1, j-1)
                    ),
                    (i-1, j),
                    (
                        (i+1, j+1), (i, j+1), (i-1, j+1)
                    ),
                    (i+1, j)
                )
    projector.(Ref(network), Ref(vertex), neighbours)
end

node_index(peps::AbstractGibbsNetwork, node::NTuple{2, Int}) = peps.ncols * (node[1] - 1) + node[2]

node_from_index(peps::AbstractGibbsNetwork, index::Int) =
    ((index-1) ÷ peps.ncols + 1, mod1(index, peps.ncols))

function boundary(peps::PEPSNetwork{false}, node::NTuple{2, Int})
    i, j = node
    x = (-4, -2)
    vcat(
        [
            [(x, x), ((i, k), (i+1, k))] for k ∈ 1:j-1
        ]...,
            (x, x), ((i, j-1), (i, j)),
        [
            [(x, x), ((i-1, k), (i, k))] for k ∈ j:peps.ncols
        ]...
    )
end

function boundary(peps::PEPSNetwork{true}, node::NTuple{2, Int})
    i, j = node
    vcat(
        [
            [((i, k-1), (i+1, k), (i, k), (i+1, k-1)), ((i, k), (i+1, k))] for k ∈ 1:j-1
        ]...,
        (
            (i, j-1), (i+1, j)
        ),
        (
            (i, j-1), (i, j)
        ),
        (
            (i-1, j-1), (i, j)
        ),
        (
            (i-1, j), (i, j)
        ),
        [
            [((i-1, k-1), (i, k), (i-1, k), (i, k-1)), ((i-1, k), (i, k))] for k ∈ j+1:peps.ncols
        ]...
    )
end

function conditional_probability(peps::PEPSNetwork{false}, w::Vector{Int})
    i, j = node_from_index(peps, length(w)+1)
    ∂v = boundary_state(peps, w, (i, j))

    L = left_env(peps, i, ∂v[1:2*j-1])
    R = right_env(peps, i, ∂v[2*j+3 : 2*peps.ncols+2])
    A = reduced_site_tensor(peps, (i, j), ∂v[2*j], ∂v[2*j+2])

    ψ = dressed_mps(peps, i)
    M = ψ[2 * j]

    @tensor prob[σ] := L[x] * M[x, d, y] * A[r, d, σ] *
                       R[y, r] order = (x, d, r, y)

    normalize_probability(prob)
end

function conditional_probability(peps::PEPSNetwork{true}, w::Vector{Int})
    i, j = node_from_index(peps, length(w)+1)
    ∂v = boundary_state(peps, w, (i, j))

    L = left_env(peps, i, ∂v[1:2*j-2])
    R = right_env(peps, i, ∂v[2*j+3 : 2*peps.ncols+2])
    A = reduced_site_tensor(peps, (i, j), ∂v[2*j-1], ∂v[2*j], ∂v[2*j+1], ∂v[2*j+2])

    ψ = dressed_mps(peps, i)
    MX, M = ψ[2*j-1], ψ[2*j]

    @tensor prob[σ] := L[x] * MX[x, m, y] * M[y, l, z] * R[z, k] *
                        A[k, l, m, σ] order = (x, y, z, k, l, m)

    normalize_probability(prob)
end

function bond_energy(
    network::AbstractGibbsNetwork,
    u::NTuple{2, Int},
    v::NTuple{2, Int},
    σ::Int
)
    fg_u, fg_v = network.vertex_map(u), network.vertex_map(v)
    if has_edge(network.factor_graph, fg_u, fg_v)
        pu, en, pv = get_prop.(Ref(network.factor_graph), Ref(fg_u), Ref(fg_v), (:pl, :en, :pr))
        energies = (pu * (en * pv[:, σ:σ]))'
    elseif has_edge(network.factor_graph, fg_v, fg_u)
        pv, en, pu = get_prop.(Ref(network.factor_graph), Ref(fg_v), Ref(fg_u), (:pl, :en, :pr))
        energies = (pv[σ:σ, :] * en) * pu
    else
        energies = zeros(length(local_energy(network, u)))
    end
    vec(energies)
end

function update_energy(network::PEPSNetwork{false}, σ::Vector{Int})
    i, j = node_from_index(network, length(σ)+1)
    bond_energy(network, (i, j), (i, j-1), local_state_for_node(network, σ, (i, j-1))) +
    bond_energy(network, (i, j), (i-1, j), local_state_for_node(network, σ, (i-1, j))) +
    local_energy(network, (i, j))
end

function update_energy(network::PEPSNetwork{true}, σ::Vector{Int})
    i, j = node_from_index(network, length(σ)+1)
    bond_energy(network, (i, j), (i, j-1), local_state_for_node(network, σ, (i, j-1))) +
    bond_energy(network, (i, j), (i-1, j), local_state_for_node(network, σ, (i-1, j))) +
    bond_energy(network, (i, j), (i-1, j-1), local_state_for_node(network, σ, (i-1, j-1))) +
    bond_energy(network, (i, j), (i-1, j+1), local_state_for_node(network, σ, (i-1, j+1))) +
    local_energy(network, (i, j))
end


text:
function test_pgbf()
  s = pgbf(1.0)
  px = pgbf(1.0,0,0,0,1,0,0)
  @assert isapprox(amplitude(s,0,0,0),0.71270547)
  @assert isapprox(amplitude(px,0,0,0),0)
  @printf "test_pgbf is passed\n"
end

function test_cgbf()
  c = cgbf(0.0,0.0,0.0)
  push!(c,1,1)
  @assert isapprox(amplitude(c,0,0,0),0.71270547)
  c2 = cgbf(0,0,0)
  push!(c2,1,0.2)
  push!(c2,0.5,0.2)
  @assert isapprox(overlap(c2,c2),1)
  @printf "test_cgbf is passed\n"
end

function test_overlap()
  s = pgbf(1.0)
  px = pgbf(1.0,0,0,0,1,0,0)
  @assert overlap1d(0,0,0.,0.,1.) == 1
  @assert gaussian_product_center(s,s) == [0,0,0]
  @assert isapprox(overlap(s,s),1)
  @assert isapprox(overlap(px,px),1)
  @assert isapprox(overlap(s,px),0)
  @assert binomial_prefactor(0,0,0,0.,0.) == 1
end

function test_utils()
  @assert factorial2(6)==48
  @assert collect(pairs(3)) == Any[(1,1),(1,2),(1,3),(2,2),(2,3),(3,3)]
  @assert collect(pairs(3,"subdiag")) == Any[(1,2),(1,3),(2,3)]
  @assert collect(pairs(2,"rect")) == Any[(1,1),(1,2),(2,1),(2,2)]
  @assert iindex(1,1,1,1) == 1
  @assert iindex(1,1,1,2) == iindex(1,1,2,1) == iindex(1,2,1,1) == iindex(2,1,1,1) == 2
  @assert iindex(1,1,2,2) == iindex(2,2,1,1) == 4
  @printf "test_utils is passed\n"
end


text:
function test_multinomial_logistic_loss_layer(backend::Backend, tensor_dim, class_weights, T, eps)
  println("-- Testing MultinomialLogisticLossLayer{$(class_weights[1]),$(class_weights[2])} on $(typeof(backend)){$T}...")

  dims = abs(rand(Int,tensor_dim)) % 6 + 6
  if class_weights[1] != :no
    op_dim = tensor_dim-1
  else
    op_dim = max(abs(rand(Int)) % tensor_dim, 1)
  end
  println("    > $dims (operate on dimension $op_dim)")

  dims_label = copy(dims); dims_label[op_dim] = 1; dims_label = tuple(dims_label...)
  dims = tuple(dims...)
  channels = dims[op_dim]

  prob = abs(rand(T, dims)) + 0.01

  label = abs(rand(Int, dims_label)) % channels
  label = convert(Array{T}, label)

  prob_blob = make_blob(backend, prob)
  label_blob = make_blob(backend, label)
  inputs = Blob[prob_blob, label_blob]

  if class_weights[1] == :equal
    weights = ones(T, channels)
  elseif class_weights[1] == :local
    weights = rand(T, channels)
  elseif class_weights[1] == :global
    weights = round(1000*rand(T, dims[1:end-1]))/1000
  else
    @assert class_weights[1] == :no
    weights = []
  end

  layer = MultinomialLogisticLossLayer(bottoms=[:pred, :labels],
      weights=weights, normalize=class_weights[2], dim=op_dim)
  state = setup(backend, layer, inputs, Blob[])

  forward(backend, state, inputs)

  if class_weights[1] == :local || class_weights[1] == :equal
    new_shape = ones(Int, tensor_dim-1); new_shape[op_dim] = dims[op_dim]
    rep_shape = [dims[1:end-1]...]; rep_shape[op_dim] = 1
    weights = repeat(reshape(weights, new_shape...), inner=rep_shape)
  end
  if class_weights[2] == :local
    weights = weights .* (channels ./ sum(weights,op_dim))
  elseif class_weights[2] == :global
    weights = weights * (prod(dims[1:end-1]) / sum(weights))
  else
    @assert class_weights[2] == :no
  end

  new_shape = [size(weights)..., 1]
  rep_shape = ones(Int, tensor_dim); rep_shape[end] = dims[end]
  weights = repeat(reshape(weights, new_shape...), inner=rep_shape)

  expected_loss = convert(T, 0)
  dim_pre, dim_prob, dim_post = split_dims(prob, op_dim)
  prob = reshape(prob, (dim_pre, dim_prob, dim_post))
  if !isempty(weights)
    weights = reshape(weights, (dim_pre, dim_prob, dim_post))
  end
  label = reshape(label, (dim_pre, 1, dim_post))
  for i = 1:dim_pre
    for j = 1:dim_post
      if isempty(weights)
        expected_loss += -log(prob[i, round(Int, label[i,1,j])+1, j])
      else
        idx = round(Int, label[i,1,j])+1
        expected_loss += -log(prob[i,idx,j]) * weights[i,idx,j]
      end
    end
  end
  expected_loss /= prod(dims) / dims[op_dim]

  @test -eps < state.loss - expected_loss < eps

  shutdown(backend, state)
end
function test_multinomial_logistic_loss_layer(backend::Backend, class_weights, T, eps)
  for i in [2,4,5]
    test_multinomial_logistic_loss_layer(backend, i, class_weights, T, eps)
  end
end
function test_multinomial_logistic_loss_layer(backend::Backend, T, eps)
  for class_weights in ((:equal,:local),(:local,:local),(:global,:global),(:global,:local),(:no,:no))
    test_multinomial_logistic_loss_layer(backend, class_weights, T, eps)
  end
end
function test_multinomial_logistic_loss_layer(backend::Backend)
  test_multinomial_logistic_loss_layer(backend, Float64, 1e-5)
  test_multinomial_logistic_loss_layer(backend, Float32, 1e-2)
end

if test_cpu
  test_multinomial_logistic_loss_layer(backend_cpu)
end
if test_gpu
  test_multinomial_logistic_loss_layer(backend_gpu)
end


text:
using RollingFunctions

data = parse.(Int, split(read(joinpath(@__DIR__, "./1.txt"), String), "\n"))

sum(map((x, y) -> x > y, data[2:end], data[1:end-1]))

sw = rolling((x)->sum(x),data,3)

sum(map((x, y) -> x > y, sw[2:end], sw[1:end-1]))


text:
abstract type AbstractCamParams end

struct RgbdCamParams <: AbstractCamParams
    width::Int64
    height::Int64
    fx::Float64
    fy::Float64
    cx::Float64
    cy::Float64

    # virtual stereo baseline
    baseline::Float64
    depth_factor::Float64

    k1::Float64
    k2::Float64
    p1::Float64
    p2::Float64
    k3::Float64

    # for calculation
    invfx::Float64
    invfy::Float64
    # Stereo baseline multiplied by fx.
    bfx::Float64
end

function RgbdCamParams(;
    width::Int64,
    height::Int64,
    fx::Float64,
    fy::Float64,
    cx::Float64,
    cy::Float64,
    depth_factor::Float64,
    baseline::Float64=0.1,
    k1::Float64=0.0,
    k2::Float64=0.0,
    p1::Float64=0.0,
    p2::Float64=0.0,
    k3::Float64=0.0,
)
    invfx = 1/fx
    invfy = 1/fy

    bfx = baseline * fx

    cam = RgbdCamParams(
        width,
        height,
        fx,
        fy,
        cx,
        cy,
        baseline,
        depth_factor,
        k1,
        k2,
        p1,
        p2,
        k3,
        invfx,
        invfy,
        bfx,
    )
    return cam
end

# stereo camera

struct StereoCamParams <: AbstractCamParams
    width::Int64
    height::Int64
    fx::Float64
    fy::Float64
    cx::Float64
    cy::Float64

    # baseline in meter
    baseline::Float64

    k1::Float64
    k2::Float64
    p1::Float64
    p2::Float64
    k3::Float64

    # for calculation
    invfx::Float64
    invfy::Float64
    # Stereo baseline multiplied by fx.
    bfx::Float64
end

function StereoCamParams(;
    width::Int64,
    height::Int64,
    fx::Float64,
    fy::Float64,
    cx::Float64,
    cy::Float64,
    baseline::Float64=0.1,
    k1::Float64=0.0,
    k2::Float64=0.0,
    p1::Float64=0.0,
    p2::Float64=0.0,
    k3::Float64=0.0,
)
    invfx = 1/fx
    invfy = 1/fy

    bfx = baseline * fx

    cam = StereoCamParams(
        width,
        height,
        fx,
        fy,
        cx,
        cy,
        baseline,
        k1,
        k2,
        p1,
        p2,
        k3,
        invfx,
        invfy,
        bfx,
    )
    return cam
end


text:
scores = [0.0, 1.0, 1.0, 0.0, 1.0, 0.0]
events = [0,1,1,1,0,0]

@test auc(scores, events) == 2/3

scores = [1.0,0.5,0.0,0.25,0.75] .- 1.0

@test all(get_quantile_scores(scores) .== scores .+ 1.0)

x = zeros(10, 3)
y = ones(10, 3)
z = fill(2.0, 10, 3)

@test all(compute_ensemble(x,z, ensemble = "mean") .== ones(10,3))
@test all(compute_ensemble(x,y, ensemble = "max") .== ones(10,3))
@test all(compute_ensemble(x,y, ensemble = "min") .== zeros(10,3))

@test all(compute_ensemble(x,y,z, ensemble = "mean") .== ones(10,3))
@test all(compute_ensemble(x,y,z, ensemble = "max") .== fill(2.0, 10, 3))
@test all(compute_ensemble(x,y,z, ensemble = "min") .== zeros(10,3))
@test all(compute_ensemble(x,y,z, ensemble = "median") .== ones(10,3))

@test all(compute_ensemble(x,y,y,z, ensemble = "mean") .== ones(10,3))
@test all(compute_ensemble(x,y,y,z, ensemble = "max") .== fill(2.0, 10, 3))
@test all(compute_ensemble(x,y,y,z, ensemble = "min") .== zeros(10,3))
@test all(compute_ensemble(x,y,y,z, ensemble = "median") .== ones(10,3))


text:
using Interpolations

export interpolatable_field

"""
    interpolatable_field(x,y,f::ScalarGridData)

Generates an interpolatable version of grid data `f`, based on coordinates in
`x` and `y` (which should be in range form). The output can be called as a function
with coordinate pairs as arguments.
"""
function interpolatable_field(x::AbstractArray,y::AbstractArray,f::ScalarGridData)

  (length(x) == size(f,1) && length(y) == size(f,2)) ||
    error("incompatible sizes of field and coordinate data")

  return CubicSplineInterpolation((x, y),f, extrapolation_bc = (Flat(),Flat()))

end

"""
    interpolatable_field(xu,yu,xv,yv,f::EdgeData/NodePair)

Generates an interpolatable version of grid data `f`, based on coordinates in
`xu`, `yu`, `xv`, `yv` (which should be in range form). The output can be called as a function
with coordinate pairs as arguments.
"""
function interpolatable_field(xu::AbstractRange,yu::AbstractRange,
                              xv::AbstractRange,yv::AbstractRange,q::Union{Edges,NodePair})

  return interpolatable_field(xu,yu,q.u), interpolatable_field(xv,yv,q.v)

end

"""
    interpolatable_field(f::GridData,g::PhysicalGrid)

Generates an interpolatable version of grid data `f`, based on grid `g`. The output can be called as a function
with coordinate pairs as arguments.
"""
interpolatable_field(f::GridData,g::PhysicalGrid) = interpolatable_field(coordinates(f,g)...,f)

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids