text:
using Logging

# from https://discourse.julialang.org/t/task-switch-not-allowed-from-inside-staged-nor-pure-functions/20488/9
# define safe loggers for use in generated functions (where task switches are not allowed)
for level in [:debug, :info, :warn, :error]
    @eval begin
        macro $(Symbol("safe_$level"))(ex...)
            macrocall = :(@placeholder $(ex...))
            # NOTE: `@placeholder` in order to avoid hard-coding @__LINE__ etc
            macrocall.args[1] = Symbol($"@$level")
            quote
                ret = nothing
                with_logger(Logging.ConsoleLogger(Base.stderr, global_logger().min_level)) do
                    ret = $(esc(macrocall))
                end
                ret
            end
        end
    end
end

const message_severities = [
    DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT,
]

"""
Default callback for debugging with [`DebugUtilsMessengerEXT`](@ref).
"""
function default_debug_callback(message_severity, message_type, callback_data_ptr, user_data_ptr)
    callback_data_ptr == C_NULL && return UInt32(0)
    callback_data = unsafe_load(callback_data_ptr)
    message = unsafe_string(callback_data.pMessage)

    # ignore messages about available device extensions
    if !startswith(message, "Device Extension: VK")
        id_name = unsafe_string(callback_data.pMessageIdName)
        msg_type = @match message_type begin
            &DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT => "General"
            &DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT => "Validation"
            &DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT => "Performance"
            _ => error("Unknown message type $message_type")
        end
        log = string("$msg_type ($id_name): $message")
        @match message_severity begin
            &DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT => @safe_debug log
            &DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT => @safe_info log
            &DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT => @safe_warn log
            &DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT => @safe_error log
            _ => error("Unknown message severity $message_severity")
        end
    end
    UInt32(0)
end

"""
Register a user-defined callback and return the corresponding messenger.
All the levels from `min_severity` will be included. Note that this controls only what messages are sent to the callback. The logging function may use logging macros such as `@info` or `@error` to easily filter logs through the Julia logging system.

A default function [`default_debug_callback`](@ref) can be converted to a function pointer to use as a callback.

!!! warning
    `callback` must be a function pointer of type `Ptr{Nothing}` obtained from a `callback_f` function as follows:
    `callback = @cfunction(callback_f, UInt32, (DebugUtilsMessageSeverityFlagBitsEXT, DebugUtilsMessageTypeFlagBitsEXT, Ptr{vk.VkDebugUtilsMessengerCallbackDataEXT}, Ptr{Cvoid}))`
    with `callback_f` a Julia function with a signature matching the `@cfunction` call.
"""
function DebugUtilsMessengerEXT(
    instance::Instance,
    callback::Ptr{Nothing};
    min_severity = DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT,
    types = DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT | DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT | DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT
    )
    index = findfirst(==(min_severity), message_severities)
    severity_bits = message_severities[index:end]
    DebugUtilsMessengerEXT(instance, |(severity_bits...), types, callback)
end


text:
# Hello world example, similar to the Boost.Python hello world

using CxxWrap
using Base.Test

# Wrap the functions defined in C++
wrap_modules(CxxWrap._l_inheritance)

using CppInheritance.A, CppInheritance.B

# Default constructor
b = B()
@test CppInheritance.message(b) == "B"

ext: jl


text:
module TestB
    c = 4
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

# misc util tests
x = sprand(m,n,0.2)
d_x = CudaSparseMatrixCSC(x)
@test device(x)   == -1
@test length(d_x) == m*n
@test size(d_x)   == (m,n)
@test size(d_x,1) == m
@test size(d_x,2) == n
y = sprand(k,n,0.2)
d_y = CudaSparseMatrixCSC(y)
@test_throws ArgumentError copy!(d_y,d_x)
d_y = CUSPARSE.switch2csr(d_y)
d_x = CUSPARSE.switch2csr(d_x)
@test_throws ArgumentError copy!(d_y,d_x)
d_y = CUSPARSE.switch2bsr(d_y,convert(Cint,blockdim))
d_x = CUSPARSE.switch2bsr(d_x,convert(Cint,blockdim))
@test_throws ArgumentError copy!(d_y,d_x)

# misc char tests

@test_throws ArgumentError CUSPARSE.cusparseop('Z')
@test_throws ArgumentError CUSPARSE.cusparsetype('Z')
@test_throws ArgumentError CUSPARSE.cusparsefill('Z')
@test_throws ArgumentError CUSPARSE.cusparsediag('Z')
@test_throws ArgumentError CUSPARSE.cusparsedir('Z')
@test_throws ArgumentError CUSPARSE.cusparseindex('A')

@test CUSPARSE.cusparseop('N') == CUSPARSE.CUSPARSE_OPERATION_NON_TRANSPOSE
@test CUSPARSE.cusparseop('C') == CUSPARSE.CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE
@test CUSPARSE.cusparseop('T') == CUSPARSE.CUSPARSE_OPERATION_TRANSPOSE
@test CUSPARSE.cusparsefill('U') == CUSPARSE.CUSPARSE_FILL_MODE_UPPER
@test CUSPARSE.cusparsefill('L') == CUSPARSE.CUSPARSE_FILL_MODE_LOWER
@test CUSPARSE.cusparsediag('U') == CUSPARSE.CUSPARSE_DIAG_TYPE_UNIT
@test CUSPARSE.cusparsediag('N') == CUSPARSE.CUSPARSE_DIAG_TYPE_NON_UNIT
@test CUSPARSE.cusparseindex('Z') == CUSPARSE.CUSPARSE_INDEX_BASE_ZERO
@test CUSPARSE.cusparseindex('O') == CUSPARSE.CUSPARSE_INDEX_BASE_ONE

# conversion
function test_make_csc(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSC(x)
    h_x = to_host(d_x)
    @test h_x == x
    @test eltype(d_x) == elty
end
test_make_csc(Float32)
test_make_csc(Float64)
test_make_csc(Complex64)
test_make_csc(Complex128)

function test_make_csr(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSR(x)
    h_x = to_host(d_x)
    @test h_x == x
end
test_make_csr(Float32)
test_make_csr(Float64)
test_make_csr(Complex64)
test_make_csr(Complex128)

function test_convert_r2c(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSR(x)
    d_x = CUSPARSE.switch2csc(d_x)
    h_x = to_host(d_x)
    @test h_x.rowval == x.rowval
    @test_approx_eq(h_x.nzval,x.nzval)
end
test_convert_r2c(Float32)
test_convert_r2c(Float64)
test_convert_r2c(Complex64)
test_convert_r2c(Complex128)

function test_convert_r2b(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSR(x)
    d_x = CUSPARSE.switch2bsr(d_x,convert(Cint,blockdim))
    d_x = CUSPARSE.switch2csr(d_x)
    h_x = to_host(d_x)
    @test_approx_eq(h_x,x)
end
test_convert_r2b(Float32)
test_convert_r2b(Float64)
test_convert_r2b(Complex64)
test_convert_r2b(Complex128)

function test_convert_c2b(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSC(x)
    d_x = CUSPARSE.switch2bsr(d_x,convert(Cint,blockdim))
    d_x = CUSPARSE.switch2csc(d_x)
    h_x = to_host(d_x)
    @test_approx_eq(h_x,x)
end
test_convert_c2b(Float32)
test_convert_c2b(Float64)
test_convert_c2b(Complex64)
test_convert_c2b(Complex128)

function test_convert_c2h(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSC(x)
    d_x = CUSPARSE.switch2hyb(d_x)
    d_y = CUSPARSE.switch2csc(d_x)
    CUSPARSE.cusparseDestroyHybMat(d_x.Mat)
    h_x = to_host(d_y)
    @test h_x.rowval == x.rowval
    @test_approx_eq(h_x.nzval,x.nzval)
end
test_convert_c2h(Float32)
test_convert_c2h(Float64)
test_convert_c2h(Complex64)
test_convert_c2h(Complex128)

function test_convert_r2h(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSR(x)
    d_x = CUSPARSE.switch2hyb(d_x)
    d_y = CUSPARSE.switch2csr(d_x)
    CUSPARSE.cusparseDestroyHybMat(d_x.Mat)
    h_x = to_host(d_y)
    @test h_x.rowval == x.rowval
    @test_approx_eq(h_x.nzval,x.nzval)
end
test_convert_r2h(Float32)
test_convert_r2h(Float64)
test_convert_r2h(Complex64)
test_convert_r2h(Complex128)

function test_convert_d2h(elty)
    x = rand(elty,m,n)
    d_x = CudaArray(x)
    d_x = CUSPARSE.sparse(d_x,'H')
    d_y = CUSPARSE.full(d_x)
    CUSPARSE.cusparseDestroyHybMat(d_x.Mat)
    h_x = to_host(d_y)
    @test_approx_eq(h_x,x)
end
test_convert_d2h(Float32)
test_convert_d2h(Float64)
test_convert_d2h(Complex64)
test_convert_d2h(Complex128)

function test_convert_d2b(elty)
    x = rand(elty,m,n)
    d_x = CudaArray(x)
    d_x = CUSPARSE.sparse(d_x,'B')
    d_y = CUSPARSE.full(d_x)
    h_x = to_host(d_y)
    @test_approx_eq(h_x,x)
end
test_convert_d2b(Float32)
test_convert_d2b(Float64)
test_convert_d2b(Complex64)
test_convert_d2b(Complex128)

function test_convert_c2r(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSC(x)
    d_x = CUSPARSE.switch2csr(d_x)
    h_x = to_host(d_x)
    @test h_x.rowval == x.rowval
    @test_approx_eq(h_x.nzval,x.nzval)
end
test_convert_c2r(Float32)
test_convert_c2r(Float64)
test_convert_c2r(Complex64)
test_convert_c2r(Complex128)

function test_convert_r2d(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSR(x)
    d_x = CUSPARSE.full(d_x)
    h_x = to_host(d_x)
    @test_approx_eq(h_x,full(x))
end
test_convert_r2d(Float32)
test_convert_r2d(Float64)
test_convert_r2d(Complex64)
test_convert_r2d(Complex128)

function test_convert_c2d(elty)
    x = sparse(rand(elty,m,n))
    d_x = CudaSparseMatrixCSC(x)
    d_x = CUSPARSE.full(d_x)
    h_x = to_host(d_x)
    @test_approx_eq(h_x,full(x))
end
test_convert_c2d(Float32)
test_convert_c2d(Float64)
test_convert_c2d(Complex64)
test_convert_c2d(Complex128)

function test_convert_d2c(elty)
    x = rand(elty,m,n)
    d_x = CudaArray(x)
    d_x = CUSPARSE.sparse(d_x,'C')
    h_x = to_host(d_x)
    @test_approx_eq(h_x,sparse(x))
end
test_convert_d2c(Float32)
test_convert_d2c(Float64)
test_convert_d2c(Complex64)
test_convert_d2c(Complex128)

function test_convert_d2r(elty)
    x = rand(elty,m,n)
    d_x = CudaArray(x)
    d_x = CUSPARSE.sparse(d_x)
    h_x = to_host(d_x)
    @test_approx_eq(h_x,sparse(x))
end
test_convert_d2r(Float32)
test_convert_d2r(Float64)
test_convert_d2r(Complex64)
test_convert_d2r(Complex128)

testnames = ["axpyi","dot","gthr","roti",
             "sctr","bsrsv","hybsv",
             "mv","mm","cssm",
             "bsrsm","gemm","geam","csic",
             "cssv","csilu","bsric",
             "bsrilu","gtsv"]
chosentests = testnames
if( !isempty(ARGS) )
    chosentests = ARGS
end

for test in chosentests
    include("$test.jl")
end


text:
###############################################################
# Part of BDisposal.jl package
###############################################################

"""
    prodIndex(gI,gO,bO,bI;retToScale,prodStructure,convexAssumption,startθ,startμ,startλ)

Compute productivity indexes

Given a set of measures of inputs, "good" ("desiderable") and "bad" ("undesiderable") outputs for different decision making units, compute their
productivity indexes improvements (or declines) between consecutive time periods.

## Parameters:
- Positional
  *  `gI`: "Good" inputs (3D matrix by DMUs, input items and periods)
  *  `gO`: "Good" outputs (3D matrix by DMUs, input items and periods)
  *  `bO`: "Bad" outputs (3D matrix by DMUs, input items and periods)
  *  `bI`: "Bad" inputs (optional 3D matrix by DMUs, input items and periods) [default: empty array]
- Keyword
  *  `retToScale`: Wheter the return to scales should be assumed "constant" (default) or "variable"
  *  `prodStructure`: Wheter the production structure should be assumed "additive" (default) or "multiplicative"
  *  `convexAssumption`: Wheter a convex production frontier should be assumed [default: `true`]

## Returns:
* A named touple where each element is a matrix of production indexes by DMUs and period passages (e.g. "year2 on year1" and "year3 on year2") or one of their decompositions.
* Currently the value reported are:
    - `prodIndexes`:          Overall production indexes
    - `prodIndexes_G`:        Decomposition for "good" inputs and outputs
    - `prodIndexes_B`:        Decomposition for "bad" inputs and outputs
    - `prodIndexes_T`:        Decomposition for the technological component, overall
    - `prodIndexes_T_O`:      Decomposition for the technological component, outputs
    - `prodIndexes_T_G_O`:    Decomposition for the technological component, good outputs
    - `prodIndexes_T_B_O`:    Decomposition for the technological component, bad outputs
    - `prodIndexes_T_I`:      Decomposition for the technological component, inputs
    - `prodIndexes_T_G_I`:    Decomposition for the technological component, good inputs
    - `prodIndexes_T_B_I`:    Decomposition for the technological component, bad inputs
    - `prodIndexes_E`:        Decomposition for the efficiency component, overall
    - `prodIndexes_E_O`:      Decomposition for the efficiency component, outputs
    - `prodIndexes_E_G_O`:    Decomposition for the efficiency component, good output
    - `prodIndexes_E_B_O`:    Decomposition for the efficiency component, bad outputs
    - `prodIndexes_E_I`:      Decomposition for the efficiency component, inputs
    - `prodIndexes_E_G_I`:    Decomposition for the efficiency component, good inputs
    - `prodIndexes_E_B_I`:    Decomposition for the efficiency component, bad inputs
    - `prodIndexes_S`:        Decomposition for the scale (residual) component, overall
    - `prodIndexes_S_O`:      Decomposition for the scale (residual) component, outputs
    - `prodIndexes_S_G_O`:    Decomposition for the scale (residual) component, good output
    - `prodIndexes_S_B_O`:    Decomposition for the scale (residual) component, bad outputs
    - `prodIndexes_S_I`:      Decomposition for the scale (residual) component, inputs
    - `prodIndexes_S_G_I`:    Decomposition for the scale (residual) component, good inputs
    - `prodIndexes_S_B_I`:    Decomposition for the scale (residual) component, bad inputs

The second and third element of the tuple are respectively the "good inputs/outputs" and "bad/inputs/outputs" components.
The first matrix can be retrieved from the two components by multiplying them (for `multiplicative` production structure)
or summing them (for `additive` production strucure).

## Description of the function

`prodIndex()` displays environmental productivity indices. These
productivity measures are implemented for different time periods (eg., years, months
etc.) or spatial units (eg., countries, cities etc.), based on the  environmental efficiency
indicators described in `efficiencyScores`. Hence, the environmental productivity indices inherit the
structure of additive and multiplicative productivity measures.

## Interpretation of the results

The additive productivity indicator shows combined desirable and undesirable outputs
productivity improvement (respectively decline) when it takes positive (respectively negative) values.
In the multiplicative context, if the productivity measure is greater (respectively lesser) than
1 then, desirable and undesirable outputs productivity increase (respectively decrease)
arises. The BDisposal package underscores the prominent sources of environmental productivity
change. The main drivers of productivity variation are technological change,
technical efficiency variation and scale efficiency change. For the additive background,
when the technological change is greater (respectively, lesser) than 0 then, technolog-
ical improvement (respectively, deterioration) occurs. A similar reasoning applies for
the additive technical and scale efficiency components. In the multiplicative context, if
the technological change is greater (respectively, lesser) than 1 then, technological in-
crease (respectively, decrrease) arises. A similar reasonnng applies for the multiplicative
technical and scale efficiency components.

## Notes:
* The decomposition by technological, efficiency and scale components is done only for the convex assumption,
as in the non-convex case the individual distance components used to compute these disaggregations are infinite.

"""
function prodIndex(gI::Array{Float64,3},gO::Array{Float64,3},bO::Array{Float64,3},bI::Array{Float64,3}=Array{Float64}(missing, (size(gI,1),0,size(gI,3)));
                   retToScale="constant",prodStructure="multiplicative",convexAssumption=true)

    nDMUs = size(gI,1)
    nPer = size(gI,3)
    nBI  = size(bI,2)

    # Output containers

    # - Global
    prodIndexes = Array{Union{Float64,Missing}}(missing, (nDMUs,nPer-1))

    # - decomposition by type of the inputs or outputs (Good vs Bad)
    prodIndexes_G = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_B = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))

    # - decomposition by tecnological/efficiency/scale (residual)
    prodIndexes_T     = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1)) # technological
    prodIndexes_T_O   = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_T_G_O = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_T_B_O = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_T_I   = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_T_G_I = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_T_B_I = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))

    prodIndexes_E     = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1)) # Efficiency
    prodIndexes_E_O   = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_E_G_O = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_E_B_O = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_E_I   = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_E_G_I = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_E_B_I = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))

    #prodIndexes_S     = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1)) # Scale (residual)
    prodIndexes_S_O   = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_S_G_O = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_S_B_O = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_S_I   = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_S_G_I = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_S_B_I = Array{Union{Float64,Missing}}(missing, (size(gI,1),nPer-1))
    prodIndexes_S     = Array{Union{Float64,Missing}}(missing, (nDMUs,nPer-1))

    noBadIndexDefault = prodStructure == "multiplicative" ? 1.0 : 0.0

    for t in 1:nPer-1
        #println(t)
        gIₜ = gI[:,:,t]
        gIᵤ = gI[:,:,t+1]
        bIₜ = bI[:,:,t]
        bIᵤ = bI[:,:,t+1]
        gOₜ = gO[:,:,t]
        gOᵤ = gO[:,:,t+1]
        bOₜ = bO[:,:,t]
        bOᵤ = bO[:,:,t+1]
        #if convexAssumption == true
            dirGI = prodStructure == "multiplicative" ? (-1,0,0,0) : (1,0,0,0)
            dirBI = prodStructure == "multiplicative" ? (0,-1,0,0) : (0,1,0,0)
            dirGO = prodStructure == "multiplicative" ? (0,0,1,0)  : (0,0,1,0)
            dirBO = prodStructure == "multiplicative" ? (0,0,0,-1) : (0,0,0,1)
        #else
        #    dirGI = prodStructure == "multiplicative" ? (1,0,0,0)  : (1,0,0,0)
        #    dirBI = prodStructure == "multiplicative" ? (0,1,0,0)  : (0,1,0,0)
        #    dirGO = prodStructure == "multiplicative" ? (0,0,-1,0) : (0,0,1,0)
        #    dirBO = prodStructure == "multiplicative" ? (0,0,0,1)  : (0,0,0,-1)
        #end
        for z in 1:nDMUs
            gIₜ₀ = gIₜ[z,:]
            bIₜ₀ = bIₜ[z,:]
            gOₜ₀ = gOₜ[z,:]
            bOₜ₀ = bOₜ[z,:]
            gIᵤ₀ = gIᵤ[z,:]
            bIᵤ₀ = bIᵤ[z,:]
            gOᵤ₀ = gOᵤ[z,:]
            bOᵤ₀ = bOᵤ[z,:]

           #if z == 1 && t == nPer-1
            #   println("here we stop and check")
           #end

            forceLinearModel = (convexAssumption == true) ?  true : false

            # t...
            idx_gi_t̃ = problem(gIᵤ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGI,crossTime=true)
            idx_gi_t = problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGI)
            idx_bi_t̃ = (nBI == 0) ? noBadIndexDefault : problem(gIₜ₀,bIᵤ₀,gOₜ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBI,crossTime=true)
            idx_bi_t = (nBI == 0) ? noBadIndexDefault : problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBI)
            idx_go_t̃ = problem(gIₜ₀,bIₜ₀,gOᵤ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGO,crossTime=true)
            idx_go_t = problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGO)
            idx_bo_t̃ = problem(gIₜ₀,bIₜ₀,gOₜ₀,bOᵤ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBO,crossTime=true)
            idx_bo_t = problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBO)

            # u...
            idx_gi_u = problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGI)
            idx_gi_ũ = problem(gIₜ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGI,crossTime=true)
            idx_bi_u = (nBI == 0) ? noBadIndexDefault : problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBI)
            idx_bi_ũ = (nBI == 0) ? noBadIndexDefault : problem(gIᵤ₀,bIₜ₀,gOᵤ₀,bOᵤ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBI,crossTime=true)
            idx_go_u = problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGO)
            idx_go_ũ = problem(gIᵤ₀,bIᵤ₀,gOₜ₀,bOᵤ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGO,crossTime=true)
            idx_bo_u = problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBO)
            idx_bo_ũ = problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOₜ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBO,crossTime=true)

            # DMU at t, frontier at t+1...
            idx_gi_tu = problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGI,crossTime=true)
            idx_bi_tu = (nBI == 0) ? noBadIndexDefault : problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBI,crossTime=true)
            idx_go_tu = problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGO,crossTime=true)
            idx_bo_tu = problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIᵤ,bIᵤ,gOᵤ,bOᵤ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBO,crossTime=true)

            # DMU at t+1, frontier at t...
            idx_gi_ut = problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGI,crossTime=true)
            idx_bi_ut = (nBI == 0) ? noBadIndexDefault : problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBI,crossTime=true)
            idx_go_ut = problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirGO,crossTime=true)
            idx_bo_ut = problem(gIᵤ₀,bIᵤ₀,gOᵤ₀,bOᵤ₀,gIₜ,bIₜ,gOₜ,bOₜ,
                        retToScale=retToScale,prodStructure=prodStructure,
                        convexAssumption=convexAssumption,forceLinearModel=forceLinearModel,
                        directions=dirBO,crossTime=true)

            # Computing aggregated for the full (gI,Bi,bO,gO) partition...
            if prodStructure == "multiplicative"
                if convexAssumption == true
                    idx_i_t = (idx_gi_t̃/idx_gi_t) * (idx_bi_t̃/idx_bi_t)
                    idx_o_t = (idx_go_t/idx_go_t̃) * (idx_bo_t/idx_bo_t̃)
                    idx_t   = idx_o_t/idx_i_t
                    idx_i_u = (idx_gi_u/idx_gi_ũ) * (idx_bi_u/idx_bi_ũ)
                    idx_o_u = (idx_go_ũ/idx_go_u) * (idx_bo_ũ/idx_bo_u)
                    idx_u   = idx_o_u/idx_i_u
                    idx     = (idx_t * idx_u)^(1/2)
                else
                    idx_i_t = (idx_gi_t/idx_gi_t̃) * (idx_bi_t/idx_bi_t̃)
                    idx_o_t = (idx_go_t/idx_go_t̃) / (idx_bo_t̃/idx_bo_t)
                    idx_t   = idx_o_t/idx_i_t
                    idx_i_u = (idx_gi_ũ/idx_gi_u) * (idx_bi_ũ/idx_bi_u)
                    idx_o_u = (idx_go_ũ/idx_go_u) / (idx_bo_u/idx_bo_ũ)
                    idx_u   = idx_o_u/idx_i_u
                    idx     = (idx_t * idx_u)^(1/2)
                end

            else # addittive
                if convexAssumption == true
                    idx_i_t = (idx_gi_t̃ - idx_gi_t) + (idx_bi_t̃ - idx_bi_t)
                    idx_o_t = (idx_go_t - idx_go_t̃) - (idx_bo_t̃ - idx_bo_t)
                    idx_t   = idx_o_t - idx_i_t
                    idx_i_u = (idx_gi_u - idx_gi_ũ) + (idx_bi_u - idx_bi_ũ)
                    idx_o_u = (idx_go_ũ - idx_go_u) - (idx_bo_u - idx_bo_ũ)
                    idx_u   = idx_o_u - idx_i_u
                    idx     = (idx_t + idx_u) / 2
                else
                    idx_i_t = -(idx_gi_t-idx_gi_t̃) - (idx_bi_t-idx_bi_t̃)
                    idx_o_t = (idx_go_t-idx_go_t̃) - (idx_bo_t̃-idx_bo_t)
                    idx_t   = idx_o_t-idx_i_t
                    idx_i_u = -(idx_gi_ũ-idx_gi_u) - (idx_bi_ũ-idx_bi_u)
                    idx_o_u = (idx_go_ũ-idx_go_u) - (idx_bo_u-idx_bo_ũ)
                    idx_u   = idx_o_u-idx_i_u
                    idx     = (idx_t + idx_u)/2
                end
            end
            prodIndexes[z,t] = idx

            if convexAssumption == true
                # Computing aggregated for the good inputs/outputs only...
                if prodStructure == "multiplicative"
                   idx_Gi_t = (idx_gi_t̃/idx_gi_t)
                   idx_Go_t = (idx_go_t/idx_go_t̃)
                   idx_Gt   = idx_Go_t/idx_Gi_t
                   idx_Gi_u = (idx_gi_u/idx_gi_ũ)
                   idx_Go_u = (idx_go_ũ/idx_go_u)
                   idx_Gu   = idx_Go_u/idx_Gi_u
                   idx_G    = (idx_Gt * idx_Gu)^(1/2)
                else
                    idx_Gi_t = (idx_gi_t̃ - idx_gi_t)
                    idx_Go_t = (idx_go_t - idx_go_t̃)
                    idx_Gt   = idx_Go_t - idx_Gi_t
                    idx_Gi_u = (idx_gi_u - idx_gi_ũ)
                    idx_Go_u = (idx_go_ũ - idx_go_u)
                    idx_Gu   = idx_Go_u - idx_Gi_u
                    idx_G     = (idx_Gt + idx_Gu) / 2
                end
                prodIndexes_G[z,t] = idx_G

                # Computing aggregated for the bad inputs/outputs only...
                if prodStructure == "multiplicative"
                   idx_Bi_t = (idx_bi_t̃/idx_bi_t)
                   idx_Bo_t = (idx_bo_t̃/idx_bo_t)
                   idx_Bt   = 1/(idx_Bo_t * idx_Bi_t)
                   idx_Bi_u = (idx_bi_u/idx_bi_ũ)
                   idx_Bo_u = (idx_bo_u/idx_bo_ũ)
                   idx_Bu   = 1/(idx_Bo_u*idx_Bi_u)
                   idx_B    = (idx_Bt * idx_Bu)^(1/2)
                else
                    idx_Bi_t = (idx_bi_t̃ - idx_bi_t)
                    idx_Bo_t = (idx_bo_t̃ - idx_bo_t)
                    idx_Bt   = - idx_Bo_t - idx_Bi_t
                    idx_Bi_u = (idx_bi_u - idx_bi_ũ)
                    idx_Bo_u = (idx_bo_u - idx_bo_ũ)
                    idx_Bu   = - idx_Bo_u - idx_Bi_u
                    idx_B    = (idx_Bt + idx_Bu) / 2
                end
                prodIndexes_B[z,t] = idx_B

                # Computing tecnological/efficiency/residual decomposition
                if prodStructure == "multiplicative"
                    # Disaggregation T/E/S...
                    idx_T_G_O = ((idx_go_t/idx_go_tu) * (idx_go_ut/idx_go_u) )^-(1/2)
                    idx_T_B_O = ((idx_bo_t/idx_bo_tu) * (idx_bo_ut/idx_bo_u) )^-(1/2)
                    idx_T_O   = idx_T_G_O * idx_T_B_O
                    idx_T_G_I = ((idx_gi_t/idx_gi_tu) * (idx_gi_ut/idx_gi_u) )^-(1/2)
                    idx_T_B_I = ((idx_bi_t/idx_bi_tu) * (idx_bi_ut/idx_bi_u) )^-(1/2)
                    idx_T_I   = idx_T_G_I * idx_T_B_I
                    idx_T     = (idx_T_O * idx_T_I)

                    idx_E_G_O = (idx_go_u/idx_go_t)^-1
                    idx_E_B_O = (idx_bo_u/idx_bo_t)^-1
                    idx_E_O   = idx_E_G_O * idx_E_B_O
                    idx_E_G_I = (idx_gi_u/idx_gi_t)^-1
                    idx_E_B_I = (idx_bi_u/idx_bi_t)^-1
                    idx_E_I   = idx_E_G_I * idx_E_B_I
                    idx_E     = (idx_E_O * idx_E_I)

                    idx_S_I   = idx / (idx_T_I * idx_E_I)
                    idx_S_O   = idx / (idx_T_O * idx_E_O)
                    idx_S_G_O = ((idx_go_t̃ / idx_go_ut) * (idx_gi_t̃ / idx_gi_t) * (idx_go_tu / idx_go_ũ) * (idx_gi_u/idx_gi_ũ)  )^-(1/2)
                    idx_S_B_O = ((idx_bo_t̃ / idx_bo_ut) * (idx_bi_t̃ / idx_bi_t) * (idx_bo_tu / idx_bo_ũ) * (idx_bi_u/idx_bi_ũ)  )^-(1/2)
                    idx_S_G_I = ((idx_gi_t̃ / idx_gi_ut) * (idx_go_t̃ / idx_go_t) * (idx_gi_tu / idx_gi_ũ) * (idx_go_u/idx_go_ũ)  )^-(1/2)
                    idx_S_B_I = ((idx_bi_t̃ / idx_bi_ut) * (idx_bo_t̃ / idx_bo_t) * (idx_bi_tu / idx_bi_ũ) * (idx_bo_u/idx_bo_ũ)  )^-(1/2)

                    idx_S     = idx / (idx_T * idx_E)

                else # additive production structure
                    # Disaggregation T/E/S...
                    idx_T_G_O = ((-idx_go_t+idx_go_tu) + (-idx_go_ut+idx_go_u) )/2
                    idx_T_B_O = ((-idx_bo_t-idx_bo_tu) + (+idx_bo_ut+idx_bo_u) )/2
                    idx_T_O   = idx_T_G_O + idx_T_B_O
                    idx_T_G_I = ((-idx_gi_t+idx_gi_tu) + (-idx_gi_ut+idx_gi_u) )/2
                    idx_T_B_I = ((-idx_bi_t+idx_bi_tu) + (-idx_bi_ut+idx_bi_u) )/2
                    idx_T_I   = idx_T_G_I + idx_T_B_I
                    idx_T     = (idx_T_O + idx_T_I)

                    idx_E_G_O = -(idx_go_u-idx_go_t)
                    idx_E_B_O = -(idx_bo_u-idx_bo_t)
                    idx_E_O   = idx_E_G_O + idx_E_B_O
                    idx_E_G_I = -(idx_gi_u-idx_gi_t)
                    idx_E_B_I = -(idx_bi_u-idx_bi_t)
                    idx_E_I   = idx_E_G_I + idx_E_B_I
                    idx_E     = (idx_E_O + idx_E_I)

                    idx_S_I   = idx - (idx_T_I + idx_E_I)
                    idx_S_O   = idx - (idx_T_O + idx_E_O)
                    idx_S_G_O = -((idx_go_t̃ - idx_go_ut) + (idx_gi_t̃ - idx_gi_t) + (idx_go_tu - idx_go_ũ) + (idx_gi_u-idx_gi_ũ)  )/2
                    idx_S_B_O = -((idx_bo_t̃ + idx_bo_ut) + (idx_bi_t̃ - idx_bi_t) + (-idx_bo_tu - idx_bo_ũ) + (idx_bi_u-idx_bi_ũ)  )/2
                    idx_S_G_I = -((idx_gi_t̃ - idx_gi_ut) + (idx_go_t̃ - idx_go_t) + (idx_gi_tu - idx_gi_ũ) + (idx_go_u-idx_go_ũ)  )/2
                    idx_S_B_I = -((idx_bi_t̃ - idx_bi_ut) + (idx_bo_t̃ - idx_bo_t) + (idx_bi_tu - idx_bi_ũ) + (idx_bo_u-idx_bo_ũ)  )/2
                    idx_S     = idx - (idx_T + idx_E)

                    #=
                    idx_T_G_O = ((idx_go_t-idx_go_ũ) + (idx_go_t̃-idx_go_u) )/2
                    idx_T_B_O = ((idx_bo_t-idx_bo_ũ) + (idx_bo_t̃-idx_bo_u) )/2
                    idx_T_O   = idx_T_G_O + idx_T_B_O
                    idx_T_G_I = ((idx_gi_t-idx_gi_ũ) + (idx_gi_t̃-idx_gi_u) )/2
                    idx_T_B_I = ((idx_bi_t-idx_bi_ũ) + (idx_bi_t̃-idx_bi_u) )/2
                    idx_T_I   = idx_T_G_I + idx_T_B_I
                    idx_T     = (idx_T_O + idx_T_I)/2

                    idx_E_G_O = idx_go_u-idx_go_t
                    idx_E_B_O = idx_bo_u-idx_bo_t
                    idx_E_O   = idx_E_G_O + idx_E_B_O
                    idx_E_G_I = idx_gi_u-idx_gi_t
                    idx_E_B_I = idx_bi_u-idx_bi_t
                    idx_E_I   = idx_E_G_I + idx_E_B_I
                    idx_E     = (idx_E_O + idx_E_I)/2

                    idx_S_G_O = idx_G -  (idx_T_G_O + idx_E_G_O)
                    idx_S_B_O = idx_B -  (idx_T_B_O + idx_E_B_O)
                    idx_S_O   = idx_S_G_O + idx_S_B_O
                    idx_S_G_I = idx_G -  (idx_T_G_I + idx_E_G_I)
                    idx_S_B_I = idx_B -  (idx_T_B_I + idx_E_B_I)
                    idx_S_I   = idx_S_G_I + idx_S_B_I
                    #idx_S     = (idx_S_O + idx_S_I)/2
                    =#
                end

            prodIndexes_T[z,t]      =  idx_T
            prodIndexes_T_O[z,t]    =  idx_T_O
            prodIndexes_T_G_O[z,t]  =  idx_T_G_O
            prodIndexes_T_B_O[z,t]  =  idx_T_B_O
            prodIndexes_T_I[z,t]    =  idx_T_I
            prodIndexes_T_G_I[z,t]  =  idx_T_G_I
            prodIndexes_T_B_I[z,t]  =  idx_T_B_I
            prodIndexes_E[z,t]      =  idx_E
            prodIndexes_E_O[z,t]    =  idx_E_O
            prodIndexes_E_G_O[z,t]  =  idx_E_G_O
            prodIndexes_E_B_O[z,t]  =  idx_E_B_O
            prodIndexes_E_I[z,t]    =  idx_E_I
            prodIndexes_E_G_I[z,t]  =  idx_E_G_I
            prodIndexes_E_B_I[z,t]  =  idx_E_B_I
            prodIndexes_S[z,t]      =  idx_S
            prodIndexes_S_O[z,t]    =  idx_S_O
            prodIndexes_S_G_O[z,t]  =  idx_S_G_O
            prodIndexes_S_B_O[z,t]  =  idx_S_B_O
            prodIndexes_S_I[z,t]    =  idx_S_I
            prodIndexes_S_G_I[z,t]  =  idx_S_G_I
            prodIndexes_S_B_I[z,t]  =  idx_S_B_I

        else # non-convex case
            # Computing aggregated for the good inputs/outputs only...
            if prodStructure == "multiplicative"
                idx_Gi_t = 1/(idx_gi_t̃*idx_gi_t)
                idx_Go_t = (idx_go_t/idx_go_t̃)
                idx_Gt   = idx_Go_t/idx_Gi_t
                idx_Gi_u = 1/(idx_gi_u/idx_gi_ũ)
                idx_Go_u = (idx_go_ũ/idx_go_u)
                idx_Gu   = idx_Go_u/idx_Gi_u
                idx_G    = (idx_Gt * idx_Gu)^(1/2)

                idx_Bi_t = 1/(idx_bi_t̃/idx_bi_t)
                idx_Bo_t = 1/(idx_bo_t̃/idx_bo_t)
                idx_Bt   = (idx_Bo_t / idx_Bi_t)
                idx_Bi_u = 1/(idx_bi_u/idx_bi_ũ)
                idx_Bo_u = 1/(idx_bo_u/idx_bo_ũ)
                idx_Bu   = (idx_Bo_u/idx_Bi_u)
                idx_B    = (idx_Bt * idx_Bu)^(1/2)

            else

                idx_Gi_t = (idx_gi_t̃-idx_gi_t)
                idx_Go_t = (idx_go_t-idx_go_t̃)
                idx_Gt   = idx_Go_t-idx_Gi_t
                idx_Gi_u = (idx_gi_u-idx_gi_ũ)
                idx_Go_u = (idx_go_ũ-idx_go_u)
                idx_Gu   = idx_Go_u-idx_Gi_u
                idx_G    = (idx_Gt + idx_Gu)/2

                idx_Bi_t = (idx_bi_t̃-idx_bi_t)
                idx_Bo_t = (idx_bo_t̃-idx_bo_t)
                idx_Bt   = (-idx_Bo_t - idx_Bi_t)
                idx_Bi_u = (idx_bi_u-idx_bi_ũ)
                idx_Bo_u = (idx_bo_u-idx_bo_ũ)
                idx_Bu   = (-idx_Bo_u-idx_Bi_u)
                idx_B    = (idx_Bt + idx_Bu)/2

            end
            prodIndexes_G[z,t] = idx_G
            prodIndexes_B[z,t] = idx_B

        end

        end # end for each DMU
    end # end for each period

    return (prodIndexes        =  prodIndexes       ,
            prodIndexes_G      =  prodIndexes_G     ,
            prodIndexes_B      =  prodIndexes_B     ,
            prodIndexes_T      =  prodIndexes_T     ,
            prodIndexes_T_O    =  prodIndexes_T_O   ,
            prodIndexes_T_G_O  =  prodIndexes_T_G_O ,
            prodIndexes_T_B_O  =  prodIndexes_T_B_O ,
            prodIndexes_T_I    =  prodIndexes_T_I   ,
            prodIndexes_T_G_I  =  prodIndexes_T_G_I ,
            prodIndexes_T_B_I  =  prodIndexes_T_B_I ,
            prodIndexes_E      =  prodIndexes_E     ,
            prodIndexes_E_O    =  prodIndexes_E_O   ,
            prodIndexes_E_G_O  =  prodIndexes_E_G_O ,
            prodIndexes_E_B_O  =  prodIndexes_E_B_O ,
            prodIndexes_E_I    =  prodIndexes_E_I   ,
            prodIndexes_E_G_I  =  prodIndexes_E_G_I ,
            prodIndexes_E_B_I  =  prodIndexes_E_B_I ,
            prodIndexes_S      =  prodIndexes_S     ,
            prodIndexes_S_O    =  prodIndexes_S_O   ,
            prodIndexes_S_G_O  =  prodIndexes_S_G_O ,
            prodIndexes_S_B_O  =  prodIndexes_S_B_O ,
            prodIndexes_S_I    =  prodIndexes_S_I   ,
            prodIndexes_S_G_I  =  prodIndexes_S_G_I ,
            prodIndexes_S_B_I  =  prodIndexes_S_B_I )
end


text:
function linsfa(signal::Array{Float64})
    w_signal = whiten!(signal)
    w_sig_diff = timediff(w_signal)
    A = w_sig_diff * w_sig_diff'
    B = w_signal * w_signal'

    # F = eigfact(A, B)
    F = eigfact(B, A)

    return F[:vectors]
end

# ensure signal is zero mean and unit variance
function whiten!(signal::Array{Float64})
    # zero mean
    signal = signal.-mean(signal,1)

    # unit variance
    signal = signal./std(signal,1)
    return signal
end

function timediff(signal::Array{Float64})
    # xdiff = diff(signal, 1)
    xdiff = [signal[1, :] ; diff(signal, 1)]
    return xdiff
end

function timediff!(signal::Array{Float64})
    signal =  [signal[1, :] ; diff(signal, 1)]
end


text:
if !isdefined(:LombScargle) using LombScargle end

function find_n_peaks_in_pgram(period::Vector, power::Vector; num_peaks::Integer = 1, exclude_period_factor::Real = 1.2 )
  @assert num_peaks>= 1
  peak_periods = zeros(num_peaks)
  peaks_found = 1
  peak_periods[peaks_found] = period[findmax(power)[2]]
  while peaks_found < num_peaks
     idx_active = trues(length(period))
     for j in 1:length(period)
     for i in 1:peaks_found
        if peak_periods[i]/exclude_period_factor <= period[j] <= peak_periods[i]*exclude_period_factor
           idx_active[j] = false
        end # if near a previous peak
     end # for over peaks
     end # for over periods
     peaks_found += 1
     peak_periods[peaks_found] = period[idx_active][findmax(power[idx_active])[2]]
  end # while more peaks to be found
  peak_periods
end

function find_n_peaks(time::Vector, obs::Vector, sigma::Vector; num_peaks::Integer = 1, exclude_period_factor::Real = 1.2 )
  pgram = lombscargle(time, obs, sigma)
  find_n_peaks_in_pgram(period(pgram), power(pgram), num_peaks=num_peaks, exclude_period_factor=exclude_period_factor )
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

import REPL
using REPL.TerminalMenus
using Test

function simulate_input(expected, menu::TerminalMenus.AbstractMenu, keys...)
    keydict =  Dict(:up => "\e[A",
                    :down => "\e[B",
                    :enter => "\r")

    for key in keys
        if isa(key, Symbol)
            write(stdin.buffer, keydict[key])
        else
            write(stdin.buffer, "$key")
        end
    end

    request(menu; suppress_output=true) == expected
end

include("radio_menu.jl")
include("multiselect_menu.jl")
include("dynamic_menu.jl")

# Legacy tests
include("legacytests/old_radio_menu.jl")
include("legacytests/old_multiselect_menu.jl")
include("legacytests/config.jl")

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
# Generic template for 1-D optimization.
# Used only to explain the parameters. This function does absolutely nothing and
# shouldn't be used. Ever.
#
# We are minimizing h(t): R → R.
#
# Inputs (common to all algorithm, specifics inputs/hyperparameters are
# explained in each function):
# h       :: Line Model.  Line model is a tool used for Line Search. We can adapt it
# 	   		 to one dimensionnal optimization. See the README to know how to create
#      		 a Line model
# nlpstop :: A Stopping object dealing with the stopping criterion of the
#    		 algorithm.
# verbose :: Bool. If true, print iteration of linesearch

function one_dimensional_template(h       :: AbstractNLPModel,
                         		  nlpstop :: AbstractStopping;
			                      verbose :: Bool = true,
								  kwargs...)
	# Core of the function

	# Outputs:
	# optimal :: Status if the problem has been solved
	# nlpstop :: Information about the output

	optimal = false;
	return optimal, nlpstop
end


text:
using FortranFiles
using Dates
using Unitful
using UnitfulAtomic

"Suffix of the file"
suffix(fname::AbstractString) = fname[findlast(isequal('.'), fname)+1:end]

"Part of the file exclude the suffix"
stem(fname::AbstractString) = fname[1:findlast(isequal('.'), fname)-1]

"Name of the file"
filename(fpath::AbstractString) = fpath[findlast(isequal('/'), fpath) + 1:end]

#const HaBohr2perAng = u"Eh_au" * u"a0_au" * u"a0_au" / u"Å" / u"ħ" * u"me_au"
#const HaBohr2perAng = u"Eh_au" * u"Å" * u"Å" / u"a0_au" / u"ħ" * u"me_au"

#Note that is an extra factor (ħ-1 me) attached to the matrix eleemnts (to make it momentum)
# For simplicity we carry it forward into the c0 term
const HaBohr2perAng = u"Eh_au" * u"a0_au" * u"a0_au" / u"Å"
const HaBohr = u"Eh_au" * u"a0_au"

"""
Read the cst_ome file with given number of bands, kpoints, and spins
"""
function read_cst_ome(fname, nb, nk, ns; endian="big-endian")
    om = Array{typeof((1.0 + im) * HaBohr2perAng)}(undef, nb, nb, 3, nk, ns)
    ffile = FortranFile(fname, convert=endian)
    for ik = 1:nk, is = 1:ns, c=1:3, ib = 1:nb, jb = 1:nb
        # Unit is Ha Bohr^2 / Ang
        @inbounds om[ib, jb, c, ik, is] = read(ffile, ComplexF64) * HaBohr2perAng
    end
    close(ffile)
    uconvert.(HaBohr, om)
end

function write_cst_ome(fname::AbstractString, om::Array{T}; endian="big-endian") where {T<:ComplexF64}
    nb, _, _, nk, ns = size(om)
    ffile = FortranFile(fname, "w", convert=endian)
    for ik = 1:nk, is = 1:ns, c=1:3, ib = 1:nb, jb = 1:nb
        # Unit is Ha Bohr^2 / Ang
        write(ffile, om[ib, jb, c, ik, is])
    end
    close(ffile)
end

function write_cst_ome(fname::AbstractString, om::Array{T}; endian="big-endian") where {T<:Quantity}
    nb, _, _, nk, ns = size(om)
    ffile = FortranFile(fname, "w", convert=endian)
    for ik = 1:nk, is = 1:ns, c=1:3, ib = 1:nb, jb = 1:nb
        # Unit is Ha Bohr^2 / Ang
        write(ffile, ustrip.(HaBohr2perAng, om[ib, jb, c, ik, is]))
    end
    close(ffile)
end

"""
    read_cst_ome(seed;)

Read the `<seed>.cst_ome` file, requires `<seed>.bands` file to be present.
"""
function read_cst_ome(seed)
    output = read_bands_castep_raw("$(seed).bands")
    ns, nk, nb = size(output.bands)
    ome = read_cst_ome("$(seed).cst_ome", nb, nk, ns)
    return ome
end

"""
    read_ome_bin(seed; legacy_format=false)

Read the `<seed>.ome_bin` file, requires `<seed>.bands` file to be present.
"""
function read_ome_bin(seed; legacy_format=false)
    output = read_bands_castep_raw("$(seed).bands")
    ns, nk, nb = size(output.bands)
    ome = read_ome_bin("$(seed).ome_bin", nb, nk, ns, legacy_format=legacy_format)
    return ome
end

"""
Read the ome_bin file
"""
function read_ome_bin(fname, nb, nk, ns; endian="big-endian", legacy_format=false)
    ffile = FortranFile(fname, convert=endian)
    version = read(ffile, Float64)
    supported = 1.0
    @assert version - supported < 0.001 "Unsupported .ome_bin file version"
    header = read(ffile, FString{80})

    # The ome_bin is written in the strange legacy format like the cst_ome file...
    if legacy_format == true
        om = Array{typeof((1.0 + im) * HaBohr2perAng)}(undef, nb, nb, 3, nk, ns)
        for ik = 1:nk, is = 1:ns, c = 1:3, ib = 1:nb, jb = 1:nb
            @inbounds om[ib, jb, c, ik, is] = read(ffile, ComplexF64) * HaBohr2perAng
        end
    else
        om = Array{typeof((1.0 + im) * HaBohr)}(undef, nb, nb, 3, nk, ns)
        for ik = 1:nk
            for is = 1:ns
                @inbounds om[:, :, :, ik, is] = read(ffile, (ComplexF64, nb, nb, 3)) * HaBohr
            end
        end
    end
    close(ffile)
    om
end

"""
    read_om_castep(ome_file, nb, nk, ns;endian="big-endian", legacy_format=false)

Read optical matrix element for CASTEP
"""
function read_om_castep(ome_file, nb, nk, ns;endian="big-endian", legacy_format=false)
    if endswith(ome_file, "cst_ome")
        out = read_cst_ome(ome_file, nb, nk, ns;endian)
    elseif endswith(ome_file, "ome_bin")
        out = read_ome_bin(ome_file, nb, nk, ns;endian, legacy_format)
    else
        throw(ErrorException("Unsupported file: $ome_file"))
    end
    out
end

"""
    read_om_castep(seed; endian="big-endian", legacy_format=false)

Read optical matrix element given a CASTEP seed name.
"""
function read_om_castep(seed; endian="big-endian", legacy_format=false)
    output = read_bands_castep_raw("$(seed).bands")
    if isfile("$(seed).ome_bin")
        ome_file = "$(seed).ome_bin"
    else
        ome_file = "$(seed).cst_ome"
    end
    ns, nk, nb = size(output.bands)
    read_om_castep(ome_file, nb, nk, ns;endian, legacy_format)
end

function write_ome_bin(fname::AbstractString, om::Array{T}; endian="big-endian", version=1.0) where {T<:ComplexF64}
    nb, _, _, nk, ns = size(om)
    ffile = FortranFile(fname, "w", convert=endian)
    write(ffile, version)
    write(ffile, FString(80, "Written by NLTools at $(now())"))
    tmp = om[:, :, :, 1, 1]
    for ik = 1:nk
        for is = 1:ns
            tmp .= om[:, :, :, ik, is]
            write(ffile, tmp)
        end
    end
    close(ffile)
end

function write_ome_bin(fname::AbstractString, om::Array{T}; endian="big-endian", version=1.0) where {T<:Quantity}
    nb, _, _, nk, ns = size(om)
    ffile = FortranFile(fname, "w", convert=endian)
    write(ffile, version)
    write(ffile, FString(80, "Written by NLTools at $(now())"))
    for ik = 1:nk
        for is = 1:ns
            # Ensure the unit written is HaBohr
            tmp = ustrip.(HaBohr, om[:, :, :, ik, is])
            write(ffile, tmp)
        end
    end
    close(ffile)
end

precompile(read_ome_bin, (String,))
precompile(read_cst_ome, (String,))


text:
using Base.Threads
using BenchmarkTools
using LinearAlgebra
using LoopVectorization
using MPI
using Strided
using UltraDark: phi_whole_step!, Grids, PencilGrids

if ~MPI.Initialized()
    MPI.Init()
end

n = 128

Δt = 1.0

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1
suite = BenchmarkGroup()

for (label, grids) in [("Grids", Grids(1.0, n)), ("PencilGrids", PencilGrids(1.0, n))]

    grids.ψx .= randn(size(grids.ψx))
    grids.ψk .= randn(size(grids.ψk))

    suite[label] = BenchmarkGroup()

    for f in [
              phi_whole_step!,
             ]
        suite[label]["$f"] = @benchmarkable $f($Δt, $grids, $nothing)
    end

end

@info "Tuning"
@show tune!(suite)

@info "Running benchmarks"
results = run(suite)

@info "Results"
@show res_min = minimum(results)
# @show res_med = median(results)

@show judge(res_min["Grids"], res_min["PencilGrids"])


text:
module TestBasics

using Test
using InitialValues
using InitialValues: isknown, hasinitialvalue

OPS = [*, +, |, &, min, max, Base.add_sum, Base.mul_prod]

@testset for op in OPS
    @test op(InitialValue(op), :anything) === :anything
    @test op(:anything, InitialValue(op)) === :anything
    @test op(INIT, :anything) === :anything
    @test op(:anything, INIT) === :anything
    @test hasinitialvalue(op)
    @test hasinitialvalue(typeof(op))
    @test isknown(InitialValue(op))
end

@testset "show" begin
    @testset "$desired" for (op, desired) in [
        (+, "InitialValue(+)"),
        (*, "InitialValue(*)"),
        (|, "InitialValue(|)"),
        (&, "InitialValue(&)"),
        (min, "InitialValue(min)"),
        (max, "InitialValue(max)"),
    ]
        @test repr(InitialValue(op); context=:limit => true) == desired
        @test repr(InitialValue(op)) == "InitialValues.$desired"
        @test string(InitialValue(op)) == "InitialValues.$desired"
    end
    @testset "INIT" begin
        @test repr(INIT; context=:limit => true) == "INIT"
        @test repr(INIT) == "InitialValues.INIT"
        @test string(INIT) == "InitialValues.INIT"
    end
end

@testset "hasinitialvalue" begin
    @test !hasinitialvalue(-)
    @test !hasinitialvalue(typeof(-))
    @test !hasinitialvalue(Int)
    @test !hasinitialvalue(Type{Int})
end

@testset "missing" begin
    @test min(InitialValue(min), missing) === missing
    @test max(InitialValue(max), missing) === missing
end

@testset "promote" begin
    for op in OPS
        T = typeof(InitialValue(op))
        @test promote_type(T, Val{0}) == Union{T,Val{0}}
    end
end

@testset "convert" begin
    @testset "float" begin
        @test float(InitialValue(+)) === 0.0
        @test float(InitialValue(*)) === 1.0
        @test float(InitialValue(Base.add_sum)) === 0.0
        @test float(InitialValue(Base.mul_prod)) === 1.0
    end
    @testset "Integer" begin
        @test Integer(InitialValue(+)) === 0
        @test Integer(InitialValue(*)) === 1
        @test Integer(InitialValue(Base.add_sum)) === 0
        @test Integer(InitialValue(Base.mul_prod)) === 1
    end
    @testset for T in [
        Int,
        UInt8,
        Float64,
        Float32,
    ]
        @test convert(T, InitialValue(+))::T == 0
        @test convert(T, InitialValue(*))::T == 1
        @test convert(T, InitialValue(Base.add_sum))::T == 0
        @test convert(T, InitialValue(Base.mul_prod))::T == 1
    end

    @test convert(String, InitialValue(*)) === ""
end

@testset "asmonoid" begin
    absmin = asmonoid() do a, b
        abs(a) < abs(b) ? a : b
    end
    @test absmin(InitialValue(absmin), Inf) === Inf
    @test absmin(missing, InitialValue(absmin)) === missing
    @test absmin(InitialValue(absmin), InitialValue(absmin)) === InitialValue(absmin)
    @test hasinitialvalue(absmin)
    @test hasinitialvalue(typeof(absmin))
    @test isknown(InitialValue(absmin))

    @test asmonoid(+) === +

    add = InitialValues.AdjoinIdentity(+)
    @test add !== (+)
    @test Base.reduce_empty(add, Int) == 0
    @test_throws Union{ArgumentError,MethodError} Base.reduce_empty(absmin, Int)
end

end  # module


text:
ways = [0, 3, 8, 9, 0]
proba = ways / sum(ways)
println(proba)


text:
using Test

@testset "Availability Exceptions" begin
    ex = AvailabilityException(1, "supply error")
    @test 1 == ex.id
    @test "supply error" == ex.message

    ex1 = AvailabilityException("EX001", "demand error")
    @test "EX001" == ex1.id
    @test "demand error" == ex1.message

end

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
using DataStructures
import Base: push!, length, keys, pop!, peek, get, delete!, haskey, empty!

mutable struct PriorityDict
    pq::PriorityQueue
    key_node_map::Dict{Any,Any}
end

PriorityDict() = PriorityDict(PriorityQueue(), Dict{Any,Any}())

length(pd::PriorityDict) = length(pd.pq)

haskey(pd::PriorityDict, key) = key in keys(pd.key_node_map)
keys(pd::PriorityDict) = keys(pd.key_node_map)

function push!(pd::PriorityDict, key, value, priority::Float64)
    if haskey(pd, key)
        error("Key already present")
    end
    entry = (key, value, priority)
    pd.key_node_map[key] = entry
    enqueue!(pd.pq, entry, priority)
    return nothing
end

function pop!(pd::PriorityDict)
    if length(pd.pq) == 0
        error("The queue is empty")
    end
    (key, value, priority) = dequeue!(pd.pq)
    delete!(pd.key_node_map, key)
    return key, value, priority
end

function peek(pd::PriorityDict)
    if length(pd.pq) == 0
        error("The queue is empty")
    end
    entry, _ = peek(pd.pq)
    return entry
end

function get(pd::PriorityDict, key, default=nothing)
    if haskey(pd, key)
        key, value, priority = pd.key_node_map[key]
        return value, priority
    else
        return default
    end
end

function delete!(pd::PriorityDict, key)
    if haskey(pd, key)
        entry = pd.key_node_map[key]
        delete!(pd.pq, entry)
        delete!(pd.key_node_map, key)
        (key, value, priority) = entry
        return value, priority
    else
        throw(KeyError(key))
    end
end

function empty!(pd::PriorityDict)
    empty!(pd.pq)
    empty!(pd.key_node_map)
end


text:
## Load the MNIST Graph and draw it
##
using DelimitedFiles
using SparseArrays
using LinearAlgebra
ei,ej = copy.(eachcol(Int.(readdlm("mnist-train-4.edges"))))
xy = readdlm("mnist-train-4.xy")
n = maximum(ei)
A = sparse(ei,ej,1,n,n)
@assert issymmetric(A)

##
using MatrixNetworks
f,lam = fiedler_vector(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=400,
  markersize=2, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.3,
  #colorbar=true,
  marker_z=f,
  axis_buffer=0.02, background=nothing)
##
savefig("mnist-train-4-fiedler.png")

##
cut = spectral_cut(A)
x = zeros(n)
x[cut.set] .= 1
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=400,
  markersize=2, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.3,
  #colorbar=true,
  marker_z=x,
  axis_buffer=0.02, background=nothing)
for i=1:100
  annotate!(xy[i,1],xy[i,2],"$(train_y[i])",color=:white)
end
##
savefig("mnist-train-4-spectral-cut.png")

## Make a figure for the label 1
# find the first 10 nodes labeled 1
using MLDatasets
train_x, train_y = MNIST.traindata()
nlabels = 10
S = findall(train_y .== 1)[1:nlabels]
## Show the labeled nodes.
scatter(xy[:,1],xy[:,2])
scatter!(xy[S,1], xy[S,2], markersize=25)
##
# Form the incidence matrix for the least squares problem
function incidence_matrix(A)
  ei,ej = findnz(triu(A,1))
  n = size(A,1)
  m = length(ei)
  return sparse(1:m,ei,1,m,n) - sparse(1:m,ej,1,m,n)
end
B = incidence_matrix(A)
IS = sparse(0.1*I, n,n)
IS[S,S] .*= 1/maximum(IS)
labels = zeros(n)
labels[S] .= 1.0
## solve the least squares problems
#x = [B; IS] \ [zeros(size(B,1)); 1]
AtA = [B; IS]'* [B; IS]
Atb = [B; IS]'* [zeros(size(B,1)); labels]
##
x = AtA\Atb
##
scatter(xy[:,1],xy[:,2], marker_z = log10.(x), markerstrokewidth=0)

## Show the labeled nodes.
scatter(xy[:,1],xy[:,2],background=nothing,framestyle=:none,markersize=2,
  size=(1200,1200),markerstrokewidth=0, alpha=0.5,dpi=300)
scatter!(xy[S,1], xy[S,2], markersize=12)
for i in S
  annotate!(xy[i,1],xy[i,2], "$(train_y[i])")
end
plot!(legend=false)
savefig("mnist-train-4-labels.png")
##
scatter(xy[:,1],xy[:,2],background=nothing,framestyle=:none,markersize=2,
  size=(1200,1200),markerstrokewidth=0, marker_z = log10.(x),
  colorbar=false, alpha=0.5,dpi=300,label="")
savefig("mnist-train-4-soln.png")


text:
abstract type LimiterFunction end
abstract type minmod <:LimiterFunction end
abstract type superbee <:LimiterFunction end
abstract type vanAlbada <:LimiterFunction end
abstract type vanAlbada2 <:LimiterFunction end
abstract type vanLeer <:LimiterFunction end
abstract type ospre <:LimiterFunction end
abstract type hquick <:LimiterFunction end
abstract type koren <:LimiterFunction end

struct FluxLimiter{T<:LimiterFunction} end

struct KnickPointsParameters{T, N}
    Kr::T
    h::T
    kappa_a::T
    m::T
    n::T
    φ::Type{FluxLimiter{N}}

    function KnickPointsParameters{T, N}(;
        Kr = 5e-7,    # whipple and tucker 1999_ dimensional coeff of erosion
        h = 1.92,    # whipple and tucker 1999
        kappa_a = 4.6071,  # whipple and tucker 1999_ area-length coeff
        m = 1.0,
        n = 1.1,
        φ = FluxLimiter{vanAlbada},
        ) where {T,N}

        new{T, N}(
            Kr,
            h,
            kappa_a,
            m,
            n,
            φ,
        )

    end

end

struct KnickPointsArrays{T}
    L::T
    r::T
    rf::T

    function KnickPointsArrays{T}(River, nn, RiverPhysics) where T
        x, z, Δx = River.x, River.x, River.Δx
        L = riverlength(x, z, nn)
        r = @. -RiverPhysics.Kr * (RiverPhysics.kappa_a^RiverPhysics.m) * (L^(RiverPhysics.h * RiverPhysics.m)) / Δx # EQUATION
        rf = zero(z)
        new{T}(L,r,rf)
    end

end

struct TerraceParameters{T}
    βz::T
    h_wb::T
    P0::T
    P_off::T

    function TerraceParameters{T}(;
        βz = 1e-6,
        h_wb = 15.0,
        P0 = 5e-5,
        P_off = 5e-2,
    ) where T
        new{T}(
            βz,
            h_wb,
            P0,
            P_off,
        )
    end
end

struct HillSlopeParameters{T}
    κa::T
    diffusion::T
    alpha_dif::T

    function HillSlopeParameters{T}(;
        κa = 5e-5,
        diffusion = 0.1 ,
        alpha_dif = 1.0 ,
    ) where T
        new{T}(
            κa,
            diffusion,
            alpha_dif,
        )
    end

end

struct SubmarineParameters{T}
    K_s::T
    λ::T

    function SubmarineParameters{T}(;
        K_s = 1e2,
        λ = 5e-4,
    ) where T
        new{T}(
            K_s,
            λ,
        )
    end

end


text:
#julia> Pkg.activate("UtilityFunctions890")
#  Activating environment at #`~/Documents/julia/Econ890/UtilityFunctions890/Project.toml`
#julia> Pkg.add("Random")
#julia> using UtilityFunctions890
# Run code above on terminal before using

using UtilityFunctions890
using Test

# @testset "UtilityFunctions890.jl" begin
using Random

# We are testing a common interface for every utility function.
function util_test(u :: AbstractUtility)
    # When you generate random numbers, always set the seed first.
    rng = MersenneTwister(23);
    @testset "$u" begin
        betaR = 1.03;
        c = 1.0 .+ rand(rng, 4,3,2);
        util = utility(u, c);
        @test size(util) == size(c);

        # inverse utility = consumption
        invUtil = inv_utility(u, util);
        @test isapprox(c, invUtil);

        # Test marginal utility by perturbation
        mu = marg_utility(u, c);
        @test all(mu .> 0.0)
        @test size(mu) == size(c)
        dc = 1e-7;
        util2 = utility(u, c .+ dc);
        mu2 = (util2 .- util) ./ dc;
        # Precision is not perfect here
        @test isapprox(mu, mu2, rtol = 1e-4);
        invMu = inv_marg_utility(u, mu);
        @test isapprox(c, invMu)

        # For consumption growth, use Euler deviation
        cGrowth = c_growth(u, betaR);
        T = 7;
        cV = cGrowth .^ (1 : T);
        devV = euler_dev(u, cV, betaR);
        @test all(abs.(devV .< 1e-5))
    end
end

@testset "Utility" begin
    for u ∈ (
        UtilityLog(),
        UtilityCRRA(2.0)
        )

        util_test(u);
    end
end


text:
module Lora

using Base.LinAlg.BLAS
using Distributions
using StatsBase
using ReverseDiffSource

import Base:
  show,
  run,
  select,
  mean,
  var,
  std
import Distributions:
  Bernoulli,
  Beta,
  Binomial,
  Cauchy,
  Exponential,
  Gamma,
  Laplace,
  LogNormal,
  Normal,
  Poisson,
  TDist,
  Uniform,
  Weibull,
  logpdf,
  logcdf,
  logccdf

export
  ### types
  MCLikelihood,
  MCChain,
  ARS,
  SliceSampler,
  MH,
  RAM,
  HMC,
  MALA,
  SMMALA,
  SerialMC,
  VanillaMCTuner,
  EmpiricalMCTuner,
  MCJob,
  PlainMCJob,
  TaskMCJob,
  ### functions
  model,
  run,
  resume,
  select,
  mean,
  mcvar,
  mcse,
  ess,
  actime,
  acceptance,
  describe,
  linearzv,
  quadraticzv

include("api/api.jl")
include("api/samples.jl")
include("api/states.jl")
include("api/chains.jl")
include("parsers/expr_funcs.jl")
include("parsers/modelparser.jl")
include("parsers/definitions/DistributionsExtensions.jl")
include("parsers/definitions/AccumulatorDerivRules.jl")
include("parsers/definitions/MCMCDerivRules.jl")
include("parsers/expr_funcs.jl")
include("models/MCLikelihood.jl")
include("models/models.jl")
include("samplers/ARS.jl")
include("samplers/SliceSampler.jl")
include("samplers/MH.jl")
include("samplers/RAM.jl")
# include("samplers/IMH.jl")
include("samplers/HMC.jl")
# include("samplers/HMCDA.jl")
# include("samplers/NUTS.jl")
include("samplers/MALA.jl")
include("samplers/SMMALA.jl")
# include("samplers/RMHMC.jl")
# include("samplers/PMALA.jl")
include("runners/SerialMC.jl")
# include("runners/SerialTempMC.jl")
# include("runners/SeqMC.jl")
include("tuners/VanillaMCTuner.jl")
include("tuners/EmpiricalMCTuner.jl")
include("jobs/PlainMCJob.jl")
include("jobs/TaskMCJob.jl")
include("jobs/jobs.jl")
include("stats/mean.jl")
include("stats/var.jl")
include("stats/zv.jl")
include("diagnostics/ess.jl")
include("diagnostics/actime.jl")
include("diagnostics/summary.jl")
# Output management: filter, merge, extract, read and write chains, convert btwn chains, arrays and dataframes
# include("output/filter.jl")
# include("output/merge.jl")
include("ui/minimal.jl")

end


text:
module Curl

export
    Easy,
        set_url,
        add_header,
        enable_progress,
        get_effective_url,
        get_response_code,
        get_response_headers,
    Multi,
        add_handle,
        remove_handle

using LibCURL
using LibCURL: curl_off_t
# not exported: https://github.com/JuliaWeb/LibCURL.jl/issues/87

include("utils.jl")

function __init__()
    @check curl_global_init(CURL_GLOBAL_ALL)
end

const CURL_VERSION = unsafe_string(curl_version())
const USER_AGENT = "$CURL_VERSION julia/$VERSION"

include("Easy.jl")
include("Multi.jl")

function add_handle(multi::Multi, easy::Easy)
    @check curl_multi_add_handle(multi.handle, easy.handle)
end

function remove_handle(multi::Multi, easy::Easy)
    @check curl_multi_remove_handle(multi.handle, easy.handle)
end

end # module


text:
"""
**Is the matrix correctly formatted?**

A *correct* matrix has only 0 and 1, two dimensions, and is square.

This function returns nothing, but raises an `AssertionError` if one of the
conditions is not met.
"""
function check_food_web(A)
  @assert size(A)[1] == size(A)[2]
  @assert length(size(A)) == 2
  @assert sum(map(x -> x ∉ [0 1], A)) == 0
end

"""
**Check initial parameters**
"""
function check_initial_parameters(parameters)
  required_keys = [
    :Z,
    :vertebrates,
    #:a_vertebrate,
    #:a_invertebrate,
    #:a_producer,
    #:m_producer,
    #:y_vertebrate,
    #:y_invertebrate,
    :e_herbivore,
    :e_carnivore,
    :h,
    #:Γ,
    :c,
    #:r,
    :K
  ]
  for k in required_keys
    @assert get(parameters, k, nothing) != nothing
  end
end

"""
**Are the simulation parameters present?**

This function will make sure that all the required parameters are here,
and that the arrays and matrices have matching dimensions.
"""
function check_parameters(parameters)
  check_initial_parameters(parameters)

  required_keys = [
    :w,
    :efficiency,
    :y,
    :x,
    #:a,
    :is_herbivore,
    :is_producer
  ]
  for k in required_keys
    @assert get(parameters, k, nothing) != nothing
  end
  @assert length(parameters[:is_producer]) == length(parameters[:is_herbivore])
  @assert size(parameters[:A]) == size(parameters[:efficiency])
  @assert length(parameters[:is_producer]) == size(parameters[:A], 1)
end

"""
TODO

"""
function check_temperature_parameters(fname::String, passed_temp_parameters::NamedTuple)
    #list the parameters required for each biological rate and for each temperature effect function
    if fname == "no_effect_x"
      required_parameters = [:a_vertebrate, :a_invertebrate, :a_producer]
    elseif fname == "no_effect_r"
      required_parameters = [:r]
    elseif fname == "no_effect_handlingt"
      required_parameters = [:y_vertebrate, :y_invertebrate]
    elseif fname == "no_effect_attackr"
      required_parameters = [:Γ]
    elseif fname == "extended_eppley_r"
      required_parameters = [:maxrate_0, :eppley_exponent, :T_opt, :range, :β]
    elseif fname == "extended_eppley_x"
      required_parameters = [:maxrate_0_producer, :maxrate_0_invertebrate, :maxrate_0_vertebrate,
                      :eppley_exponent_producer, :eppley_exponent_invertebrate, :eppley_exponent_vertebrate,
                      :T_opt_producer, :T_opt_invertebrate, :T_opt_vertebrate,
                      :range_producer, :range_invertebrate, :range_vertebrate,
                      :β_producer, :β_invertebrate, :β_vertebrate]
    elseif fname == "exponential_BA_r"
      required_parameters = [:norm_constant, :activation_energy, :T0, :β, :k, :T0K]
    elseif fname == "exponential_BA_x"
      required_parameters = [:norm_constant_producer, :norm_constant_invertebrate, :norm_constant_vertebrate,
                     :activation_energy_producer, :activation_energy_invertebrate, :activation_energy_vertebrate,
                     :T0_producer, :T0_invertebrate, :T0_vertebrate,
                     :β_producer, :β_invertebrate, :β_vertebrate]
    elseif fname == "exponential_BA_functionalr"
      required_parameters = [:norm_constant_invertebrate, :norm_constant_vertebrate,
                  :activation_energy_invertebrate, :activation_energy_vertebrate,
                  :T0_invertebrate, :T0_vertebrate,
                  :β_producer, :β_invertebrate, :β_vertebrate]
    elseif fname == "extended_BA_r"
      required_parameters = [:norm_constant, :activation_energy, :T0, :β, :deactivation_energy]
    elseif fname == "extended_BA_x"
      required_parameters = [:norm_constant_producer, :norm_constant_invertebrate, :norm_constant_vertebrate,
                     :activation_energy_producer, :activation_energy_invertebrate, :activation_energy_vertebrate,
                     :deactivation_energy_producer, :deactivation_energy_invertebrate, :deactivation_energy_vertebrate,
                     :T_opt_producer, :T_opt_invertebrate, :T_opt_vertebrate,
                     :β_producer, :β_invertebrate, :β_vertebrate]
    elseif fname == "extended_BA_attackr"
      required_parameters = [:norm_constant_invertebrate, :norm_constant_vertebrate,
                     :activation_energy_invertebrate, :activation_energy_vertebrate,
                     :deactivation_energy_invertebrate, :deactivation_energy_vertebrate,
                     :T_opt_invertebrate, :T_opt_vertebrate,
                     :β_producer, :β_invertebrate, :β_vertebrate]
    elseif fname == "gaussian_r"
      required_parameters = [:shape, :norm_constant, :range, :T_opt, :β]
    elseif fname == "gaussian_x"
      required_parameters = [:norm_constant_producer, :norm_constant_invertebrate, :norm_constant_vertebrate,
                     :T_opt_producer, :T_opt_invertebrate, :T_opt_vertebrate,
                     :β_producer, :β_invertebrate, :β_vertebrate,
                     :range_producer, :range_invertebrate, :range_vertebrate]
    elseif fname == "gaussian_functionalr"
      required_parameters = [:shape,
                  :norm_constant_invertebrate, :norm_constant_vertebrate,
                  :T_opt_invertebrate, :T_opt_vertebrate,
                  :β_producer, :β_invertebrate, :β_vertebrate]
    end
    # assert that there are no additional parameter
    for temperature_parameter in collect(keys(passed_temp_parameters)) ; @assert temperature_parameter ∈ required_parameters ; end
    # assert that each parameter has a value
    for required_temp_parameters in required_parameters ; @assert get(passed_temp_parameters, required_temp_parameters, nothing) != nothing ; end
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
export MPIAcousticPropagatorParams, MPIAcousticSource, MPIAcousticReceiver, MPIAcousticPropagator, compute_PML_Params!,
        MPIAcousticPropagatorSolver, MPISimulatedObservation!, MPIAcousticReceiver
@with_kw mutable struct MPIAcousticPropagatorParams
    # number of grids along x,y axis and time steps
    NX::Int64 = 0
    NY::Int64 = 0
    M::Int64 = 0
    N::Int64 = 0
    n::Int64 = 300
    NSTEP::Int64 = 2000 * 2

    # size of grid cell and time step
    DELTAX::Float64 = 10.
    DELTAY::Float64 = 10.
    DELTAT::Float64 = 2.e-3 / 2

    # PML boundary conditon
    USE_PML_XMIN::Bool = true
    USE_PML_XMAX::Bool = true
    USE_PML_YMIN::Bool = true
    USE_PML_YMAX::Bool = true
    NPOINTS_PML::Int64 = 12
    NPOWER::Int64 = 2
    damping_x::Union{Missing,Float64} = missing
    damping_y::Union{Missing,Float64} = missing
    Rcoef::Float64 = 0.001 # Relative reflection coefficient
    vp_ref::Float64 = 1000.

    # Auxilliary Data
    Σx::Array{Float64} = []
    Σy::Array{Float64} = []
    IJ::Array{Int64} = []
    IJn::Array{Int64}  = []
    IJp::Array{Int64} = []
    IpJ::Array{Int64} = []
    IpJp::Array{Int64} = []
    IpJn::Array{Int64} = []
    InJ::Array{Int64} = []
    InJn::Array{Int64} = []
    InJp::Array{Int64} = []

    # display params
    IT_DISPLAY::Int64 = 0

    # mpi
    II::Int64 = -1
    JJ::Int64 = -1
    PropagatorKernel::Int64 = 0 # 0 - vanilla, 1 - optimized, 2 - reference
end

"""
    MPIAcousticSource

The input is the global index.

`srci`, `srcj`: nsrc
`srcv`: nt × nsrc, each row corresponds to a source time function
"""
mutable struct MPIAcousticSource
    srci::Union{Missing, PyObject}
    srcj::Union{Missing, PyObject}
    srcv::Union{Missing, PyObject}
    function MPIAcousticSource(param::MPIAcousticPropagatorParams, srci::Array{Int64,1}, srcj::Array{Int64,1}, srcv::Union{Array{Float64,2}, PyObject})
        II, JJ = param.II, param.JJ
        nsrc = length(srci)
        srcv = constant(srcv)
        srcv = srcv'
        local_srci = Int64[]
        local_srcj = Int64[]
        local_srcv = PyObject[]
        for i = 1:nsrc
            if ((II - 1) * param.n + 1 <= srci[i] <= II * param.n) &&
                    ((JJ - 1) * param.n + 1 <= srcj[i] <= JJ * param.n)
                    push!(local_srci, srci[i] - (II - 1) * param.n)
                    push!(local_srcj, srcj[i] - (JJ - 1) * param.n)
                    push!(local_srcv, srcv[i])
            end
        end
        local_srcv = hcat(local_srcv...)
        if length(local_srci)>0
            local_srci, local_srcj, local_srcv = convert_to_tensor([local_srci, local_srcj, local_srcv], [Int64, Int64, Float64])
            new(local_srci, local_srcj, local_srcv)
        else
            new(missing, missing, missing)
        end
    end
end

mutable struct MPIAcousticReceiver
    rcvi::Union{Missing,Array{Int64, 1}}
    rcvj::Union{Missing,Array{Int64, 1}}
    rcvv::Union{Missing, PyObject}
    function MPIAcousticReceiver(param::MPIAcousticPropagatorParams, rcvi::Array{Int64,1}, rcvj::Array{Int64,1})
        II, JJ = param.II, param.JJ
        nrcv = length(rcvi)
        local_rcvi = Int64[]
        local_rcvj = Int64[]
        for i = 1:nrcv
            if ((II - 1) * param.n + 1 <= rcvi[i] <= II * param.n) &&
                    ((JJ - 1) * param.n + 1 <= rcvj[i] <= JJ * param.n)
                    push!(local_rcvi, rcvi[i] - (II - 1) * param.n)
                    push!(local_rcvj, rcvj[i] - (JJ - 1) * param.n)
            end
        end
        if length(local_rcvi)>0
            new(local_rcvi, local_rcvj, missing)
        else
            new(missing, missing, missing)
        end
    end
end

"""
    AcousticPropagator

A structure holding all the information for acoustic simulation
"""
mutable struct MPIAcousticPropagator
    param::MPIAcousticPropagatorParams
    src::MPIAcousticSource
    u::PyObject
    # auxillary variables for PML
    φ::PyObject
    ψ::PyObject
end

function pml_helper(x::Float64, y::Float64, param::MPIAcousticPropagatorParams)
    outx = 0.0; outy = 0.0
    ξx = param.damping_x
    Lx = param.NPOINTS_PML * param.DELTAX
    if x<Lx && param.USE_PML_XMIN
        d = abs(Lx-x)
        outx = ξx * (d/Lx - sin(2π*d/Lx)/(2π))
    elseif x>param.DELTAX*(param.NX+1)-Lx && param.USE_PML_XMAX
        d = abs(x-(param.DELTAX*(param.NX+1)-Lx))
        outx = ξx * (d/Lx - sin(2π*d/Lx)/(2π))
    end

    ξy = param.damping_y
    Ly = param.NPOINTS_PML * param.DELTAY
    if y<Ly && param.USE_PML_YMIN
        d = abs(Ly-y)
        outy = ξy * (d/Ly - sin(2π*d/Ly)/(2π))
    elseif y>param.DELTAY*(param.NY+1)-Ly && param.USE_PML_YMAX
        d = abs(y-(param.DELTAY*(param.NY+1)-Ly))
        outy = ξy * (d/Ly - sin(2π*d/Ly)/(2π))
    end

    return outx, outy
end

function compute_PML_Params!(param::MPIAcousticPropagatorParams; check_mpi_size::Bool = true)
    NX, NY = param.NX, param.NY
    n = param.n
    r = mpi_rank()
    param.M = div(NX,n)
    param.N = div(NY,n)
    param.II, param.JJ = div(r, param.N)+1, mod(r, param.N)+1
    if check_mpi_size
        @assert mpi_size() == param.M * param.N
    end
    @assert NX>0 && NY>0
    @assert mod(NX, n)==0 && mod(NY, n)==0

    # computing damping coefficient
    c, R = param.vp_ref, param.Rcoef
    Lx = param.NPOINTS_PML * param.DELTAX
    Ly = param.NPOINTS_PML * param.DELTAY
    param.damping_x = c/Lx*log(1/R)
    param.damping_y = c/Ly*log(1/R)
    # @show c, Lx, log(1/R), param.damping_x, param.damping_y

    X = (0:n+1)*param.DELTAX .+ n*param.DELTAX*(param.II-1)
    Y = (0:n+1)*param.DELTAY .+ n*param.DELTAY*(param.JJ-1)
    Σx = zeros(n+2, n+2)
    Σy = zeros(n+2, n+2)
    for i = 1:n+2
        for j = 1:n+2
            Σx[i,j], Σy[i,j] = pml_helper(X[i], Y[j], param)
        end
    end

    param.Σx = Σx
    param.Σy = Σy
    param.IJ = get_mpi_id(2:n+1, 2:n+1, n)
    param.IJn = get_mpi_id(2:n+1, 1:n, n)
    param.IJp = get_mpi_id(2:n+1, 3:n+2,  n)
    param.IpJ = get_mpi_id(3:n+2, 2:n+1, n)
    param.IpJp = get_mpi_id(3:n+2, 3:n+2, n)
    param.IpJn = get_mpi_id(3:n+2, 1:n, n)
    param.InJ = get_mpi_id(1:n, 2:n+1, n)
    param.InJn = get_mpi_id(1:n, 1:n, n)
    param.InJp = get_mpi_id(1:n, 3:n+2, n)
    return param
end

function get_mpi_id(a, b, n)
    idx = Int64[]
    for i = 1:length(a)
        for j = 1:length(b)
            push!(idx, (n+2)*(a[i]-1)+b[j])
        end
    end
    idx
end

function one_step(param::MPIAcousticPropagatorParams, w::PyObject, wold::PyObject,
        φ, ψ, σ::PyObject, τ::PyObject, c::PyObject, i::PyObject)
    n = param.n
    Δt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY
    IJ, IpJ, InJ, IJp, IJn, IpJp, IpJn, InJp, InJn =
        param.IJ, param.IpJ, param.InJ, param.IJp, param.IJn, param.IpJp, param.IpJn, param.InJp, param.InJn

    c = reshape(c, (-1,))
    w = reshape(mpi_halo_exchange(w, param.M, param.N, tag = 5*i, deps=cast(Float64, i)), (-1,))
    wold = reshape(mpi_halo_exchange(wold, param.M, param.N, tag = 5*i+1, deps=w[1]), (-1,))
    φ = reshape(mpi_halo_exchange(φ, param.M, param.N, tag = 5*i+2, deps=wold[1]), (-1,))
    ψ = reshape(mpi_halo_exchange(ψ, param.M, param.N, tag = 5*i+3, deps=φ[1]), (-1,))

    u = (2 - σ[IJ]*τ[IJ]*Δt^2 - 2*Δt^2/hx^2 * c - 2*Δt^2/hy^2 * c) * w[IJ] +
            c * (Δt/hx)^2  *  (w[IpJ]+w[InJ]) +
            c * (Δt/hy)^2  *  (w[IJp]+w[IJn]) +
            (Δt^2/(2hx))*(φ[IpJ]-φ[InJ]) +
            (Δt^2/(2hy))*(ψ[IJp]-ψ[IJn]) -
            (1 - (σ[IJ]+τ[IJ])*Δt/2) * wold[IJ]
    u_local = u / (1 + (σ[IJ]+τ[IJ])/2*Δt)

    u_local_nxn = reshape(u_local, (n, n))

    u = reshape(mpi_halo_exchange(u_local_nxn, param.M, param.N, tag = 5*i+4, deps=ψ[1]), (-1,))
    φ = (1. -Δt*σ[IJ]) * φ[IJ] + Δt * c * (τ[IJ] -σ[IJ])/2hx *
        (u[IpJ]-u[InJ])
    ψ = (1. -Δt*τ[IJ]) * ψ[IJ] + Δt * c * (σ[IJ] -τ[IJ])/2hy *
        (u[IJp]-u[IJn])
    φ = reshape(φ, (n, n))
    ψ = reshape(ψ, (n, n))
    u = reshape(u[IJ], (n, n))

    u, φ, ψ
end

function one_step_customop(param::MPIAcousticPropagatorParams, w::PyObject, wold::PyObject,
    φ, ψ, σ::PyObject, τ::PyObject, c::PyObject, i::PyObject)
    n = param.n
    Δt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY
    IJ, IpJ, InJ, IJp, IJn, IpJp, IpJn, InJp, InJn =
        param.IJ, param.IpJ, param.InJ, param.IJp, param.IJn, param.IpJp, param.IpJn, param.InJp, param.InJn

    c = reshape(c, (-1,))
    w = reshape(mpi_halo_exchange(w, param.M, param.N, tag = 5*i, deps=cast(Float64, i)), (-1,))
    wold = reshape(mpi_halo_exchange(wold, param.M, param.N, tag = 5*i+1, deps=w[1]), (-1,))
    φ = reshape(mpi_halo_exchange(φ, param.M, param.N, tag = 5*i+2, deps=wold[1]), (-1,))
    ψ = reshape(mpi_halo_exchange(ψ, param.M, param.N, tag = 5*i+3, deps=φ[1]), (-1,))

    phi = φ
    psi = ψ
    sigma = σ
    tau = τ
    dt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY
    nx, ny = n, n
    mpi_acoustic_one_step_ = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic","mpi_acoustic_one_step", multiple=true)
    w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny = convert_to_tensor(Any[w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny], [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64,Int64])
    u, φ, ψ = mpi_acoustic_one_step_(w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny)
    u, φ, ψ = set_shape(u, n*n), set_shape(φ, n*n), set_shape(ψ, n*n)

    # u = (2 - σ[IJ]*τ[IJ]*Δt^2 - 2*Δt^2/hx^2 * c - 2*Δt^2/hy^2 * c) * w[IJ] +
    #         c * (Δt/hx)^2  *  (w[IpJ]+w[InJ]) +
    #         c * (Δt/hy)^2  *  (w[IJp]+w[IJn]) +
    #         (Δt^2/(2hx))*(φ[IpJ]-φ[InJ]) +
    #         (Δt^2/(2hy))*(ψ[IJp]-ψ[IJn]) -
    #         (1 - (σ[IJ]+τ[IJ])*Δt/2) * wold[IJ]
    # u = u / (1 + (σ[IJ]+τ[IJ])/2*Δt)
    # φ = (1. -Δt*σ[IJ]) * φ[IJ] + Δt * c * (τ[IJ] -σ[IJ])/2hx *  (w[IpJ]-w[InJ])
    # ψ = (1. -Δt*τ[IJ]) * ψ[IJ] + Δt * c * (σ[IJ] -τ[IJ])/2hy * (w[IJp]-w[IJn])

    φ = reshape(φ, (n, n))
    ψ = reshape(ψ, (n, n))
    u = reshape(u, (n, n))
    u, φ, ψ
end

function one_step_ref(param::MPIAcousticPropagatorParams, w::PyObject, wold::PyObject,
    φ, ψ, σ::PyObject, τ::PyObject, c::PyObject, i::PyObject)
    n = param.n
    Δt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY
    IJ, IpJ, InJ, IJp, IJn, IpJp, IpJn, InJp, InJn =
        param.IJ, param.IpJ, param.InJ, param.IJp, param.IJn, param.IpJp, param.IpJn, param.InJp, param.InJn

    c = reshape(c, (-1,))
    w = reshape(mpi_halo_exchange(w, param.M, param.N, tag = 5*i, deps=cast(Float64, i)), (-1,))
    wold = reshape(mpi_halo_exchange(wold, param.M, param.N, tag = 5*i+1, deps=w[1]), (-1,))
    φ = reshape(mpi_halo_exchange(φ, param.M, param.N, tag = 5*i+2, deps=wold[1]), (-1,))
    ψ = reshape(mpi_halo_exchange(ψ, param.M, param.N, tag = 5*i+3, deps=φ[1]), (-1,))

    u = (2 - σ[IJ]*τ[IJ]*Δt^2 - 2*Δt^2/hx^2 * c - 2*Δt^2/hy^2 * c) * w[IJ] +
            c * (Δt/hx)^2  *  (w[IpJ]+w[InJ]) +
            c * (Δt/hy)^2  *  (w[IJp]+w[IJn]) +
            (Δt^2/(2hx))*(φ[IpJ]-φ[InJ]) +
            (Δt^2/(2hy))*(ψ[IJp]-ψ[IJn]) -
            (1 - (σ[IJ]+τ[IJ])*Δt/2) * wold[IJ]
    u = u / (1 + (σ[IJ]+τ[IJ])/2*Δt)
    φ = (1. -Δt*σ[IJ]) * φ[IJ] + Δt * c * (τ[IJ] -σ[IJ])/2hx *  (w[IpJ]-w[InJ])
    ψ = (1. -Δt*τ[IJ]) * ψ[IJ] + Δt * c * (σ[IJ] -τ[IJ])/2hy * (w[IJp]-w[IJn])

    φ = reshape(φ, (n, n))
    ψ = reshape(ψ, (n, n))
    u = reshape(u, (n, n))
    u, φ, ψ
end

"""
    MPIAcousticPropagatorSolver(param::MPIAcousticPropagatorParams, src::MPIAcousticSource, c::Union{PyObject, Array{Float64, 2}})

`MPIAcousticPropagatorSolver` is a MPI-enabled solver. The computational domain has a size `NX x NY` and is divided into `M x N`
"""
function MPIAcousticPropagatorSolver(param::MPIAcousticPropagatorParams, src::MPIAcousticSource, c::Union{PyObject, Array{Float64, 2}})

    c = constant(c)
    n = param.n

    @assert size(c, 1)==size(c,2)==n
    if !ismissing(src.srcv)
        @assert size(src.srcv, 2)==length(src.srci)==length(src.srcj)
        @assert size(src.srcv, 1)==param.NSTEP
    end

    σij = constant(param.Σx) # (n+2) × (n+2)
    τij = constant(param.Σy) # (n+2) × (n+2)

    σij = reshape(σij, (-1,))
    τij = reshape(τij, (-1,))

    function condition(i, ta, tφ, tψ)
        i<=param.NSTEP+1
    end

    function body(i, ta, tφ, tψ)

        if param.PropagatorKernel==0
            one_step_ = one_step
        elseif param.PropagatorKernel==1
            @info "using custom op implementation"
            one_step_ = one_step_customop
        elseif param.PropagatorKernel==2
            @info "using reference implementation"
            one_step_ = one_step_ref
        else
            error("Not implemented")
        end
        u, φ, ψ = one_step_(param, read(ta, i-1), read(ta, i-2), read(tφ, i-1), read(tψ, i-1), σij, τij, c, i)
        srci, srcj, srcv = AcousticSourceAtTimeT(src, i-1)

        if param.IT_DISPLAY>0
            op = tf.cond(tf.equal(tf.math.floormod(i,param.IT_DISPLAY),1), ()->tf.print("time step = ", i-1, " out of $(param.NSTEP)"), ()->tf.no_op())
            i = bind(i, op)
        end

        if !ismissing(srci)
            src_index = (srci - 1) * n + srcj
            u = reshape(u, (-1,))
            u = scatter_add_op(u, src_index, srcv*param.DELTAT^2)
            u = reshape(u, (n, n))
        end

        ta_, tφ_, tψ_ = write(ta, i, u), write(tφ, i, φ), write(tψ, i, ψ)
        i+1, ta_, tφ_, tψ_
    end

    tu = TensorArray(param.NSTEP+1; clear_after_read=false)
    tφ = TensorArray(param.NSTEP+1; clear_after_read=true)
    tψ = TensorArray(param.NSTEP+1; clear_after_read=true)
    tu = write(tu, 1, constant(zeros(n, n)))
    tφ = write(tφ, 1, constant(zeros(n, n)))
    tψ = write(tψ, 1, constant(zeros(n, n)))
    tu = write(tu, 2, constant(zeros(n, n)))
    tφ = write(tφ, 2, constant(zeros(n, n)))
    tψ = write(tψ, 2, constant(zeros(n, n)))
    i = constant(3, dtype=Int32)
    _, tu, tφ, tψ = while_loop(condition, body, [i,tu,tφ,tψ]; parallel_iterations=1)
    tu = set_shape(stack(tu), (param.NSTEP+1, n, n))
    tφ = set_shape(stack(tφ), (param.NSTEP+1, n, n))
    tψ = set_shape(stack(tψ), (param.NSTEP+1, n, n))

    MPIAcousticPropagator(param, src, tu, tφ, tψ)
end

function AcousticSourceAtTimeT(src::MPIAcousticSource,i::PyObject)
    if ismissing(src.srci)
        return missing, missing, missing
    else
        return src.srci, src.srcj, src.srcv[i]
    end
end

"""
    SimulatedObservation!(ap::AcousticPropagator, rcv::AcousticReceiver)

Extract and save simulated displacement u into rcv::AcousticReceiver.
"""
function MPISimulatedObservation!(ap::MPIAcousticPropagator, rcv::MPIAcousticReceiver)
    if ismissing(rcv.rcvi)
        op = group(ap.u)
        rcv.rcvv = op
        return
    end
    idx = @. (rcv.rcvi - 1) *  ap.param.n + rcv.rcvj
    @assert size(ap.u,1)==ap.param.NSTEP+1
    @assert size(ap.u,2)==ap.param.n
    rcv.rcvv = reshape(ap.u, (ap.param.NSTEP+1, -1))[:, idx]
end


text:
"""
 Module Name    : Primary solutions by heuristics functions
 Description    : The module focus on delivering a primary solution upon its
                  calling.
 Authors        : Andrew Rosemberg and Bianca Lacê

 Initial date   : 04/12/2017
"""
module feasible_solution
using JuMP

include("branch_and_bound.jl")
using Combinatorics

"""
 Function Name   : Grasp generic mip problems
 Description     : "Greedy ramdomized adaptive search procedure" porpouse is
                    to find a primary solution to a
                    Linear Integer programim problem.
                    This inplementation was optimized to solve the
                    travelling salesman problem (TSP).

"""
function _grasp_tsp(m::JuMP.Model,maxiter::Int = 200, α::Float64 = 0.25, max_iter_improv::Int= 200,max_improvs_aux::Int=200)
    num_variables = sum(m.colCat .== :Bin)
    num_nodes = Int(sqrt(num_variables))
    num_rlc = floor(num_nodes*α)
    num_rlc = Int(min(max(1,num_rlc),num_nodes))
    c = JuMP.prepAffObjective(m)
    range_nodes = 1:num_nodes
    indVariable = collect(1:m.numCols)[m.colCat .== :Bin]

    tour = fill(-1, num_nodes)

    ## first solution ##
    # restrictions
    for i = Int.(range_nodes)
        if m.objSense == :Min
            c[i + num_nodes*(i-1)] = maximum(c[1 + num_nodes*(i-1):num_nodes*(i)])  # so it wont be selected as best cadidate
        end
    end
    for i = num_variables
        if m.colUpper[i] == 0
            row = Int(floor(i/num_nodes))
            col = Int(i - num_nodes*row)
            if m.objSense == :Min
                c[row*num_nodes+col] = maximum(c[1 + num_nodes*row:num_nodes*(row+1)]) # so it wont be selected as best cadidate
            end
        end
    end

    # rlc
    # c_aux = c[1 + num_nodes*(i-1):num_nodes*(i)]
    # sortperm(c[1 + num_nodes*(i-1):num_nodes*(i)], rev=true)
    ind_aux = collect(range_nodes)
    # ind = m.objSense == :Max ? ind_aux[c_aux .>= (1-α)*maximum(c_aux)] : ind_aux[c_aux .<= (α+1)*minimum(c_aux)]

    # first iteraction
    itr = 1
    i = ind_aux[rand(1:size(ind_aux,1))]

    tour[itr] = i

    # remaining iteractions
    proh = Int.(-1*ones(num_nodes,num_nodes))
    itr+=1
    while itr <=num_nodes-1
        rlc = m.objSense == :Max ? setdiff(setdiff(sortperm(c[1 + num_nodes*(i-1):num_nodes*(i)], rev=true),tour), proh[itr,:]):
            setdiff(setdiff(sortperm(c[1 + num_nodes*(i-1):num_nodes*(i)]),tour), proh[itr,:])
        rlc = rlc[1:min(num_rlc,size(rlc,1))]
        if size(rlc,1) == 0
            itr-=1
            proh[itr,tour[itr]] = tour[itr]
            proh[itr+1,:] .=  Int.(-1*ones(num_nodes))
            tour[itr] = -1
        else
            ind = rlc[rand(1:size(rlc,1))]
            i = ind
            tour[itr] = i
            if size(setdiff(collect(1:num_nodes),tour),1) == 0
                proh[itr,tour[itr]] = tour[itr]
                proh[itr+1,:] .=  Int.(-1*ones(num_nodes))
                tour[itr] = -1
                itr-=1
            else
                itr+=1
            end
        end
    end
    i = setdiff(collect(1:num_nodes),tour)[1]
    tour[num_nodes] = i

    # check feasibility
    sol = feasible_solution._prep_solution(m,tour)
    # check feasability

    status,obval = feasible_solution._check_feasability_solution(m,sol,indVariable)
    if status == false
        # repair
    end
    ## end first solution ##
    ## ajustments solution ##
    # Initialization
    m_aux = deepcopy(m)
    max_improvs = max(num_nodes^2,max_improvs_aux)
    best_sol = deepcopy(m.colVal)
    best_obj = deepcopy(obval)
    itr_improv = 0
    itr = 1
    neighb_param1 = 2
    neighb_param2 = 1  # shift

    neighb = collect(combinations(1:num_nodes,neighb_param1))
    while itr_improv <= max_iter_improv
        # permutation
        inds = neighb[itr] #rand(1:num_nodes,2) collect(1:5)
        tour_aux = deepcopy(tour)
        tour_aux2 = deepcopy(tour[inds])
        tour_aux[inds] = circshift(tour_aux2,neighb_param2)

        # tour_aux[inds[1]] = tour[inds[2]]
        # tour_aux[inds[2]] = tour[inds[1]]

        sol = feasible_solution._prep_solution(m_aux,tour_aux)
        status,obval = feasible_solution._check_feasability_solution(m_aux,sol,indVariable)
        if status == true
            if m.objSense == :Max && obval > best_obj
                best_sol = deepcopy(m_aux.colVal)
                best_obj = deepcopy(obval)
            elseif  m.objSense == :Min && obval < best_obj
                best_sol = deepcopy(m_aux.colVal)
                best_obj = deepcopy(obval)
            end
        end
        if itr == size(neighb,1)
            if neighb_param1 == num_nodes
                neighb_param1 = 1
                neighb_param2 +=1
            end
            itr = 0
            neighb_param1 +=1
            neighb = collect(combinations(1:num_nodes,neighb_param1))
        end
        itr+=1
        itr_improv+=1
    end

    m.colVal = best_sol
    m.objVal = best_obj
    # return status
    return :SubOptimal
end

"""
Function Name   : prep solution
Description     : constructs a solution from a tour.
"""
function _prep_solution(m::JuMP.Model,tour)
    num_variables = sum(m.colCat .== :Bin)
    num_nodes = Int(sqrt(num_variables))

    # prep solution
    sol = zeros(num_variables)
    indVariable = collect(1:m.numCols)[m.colCat .== :Bin]
    row = tour[num_nodes]
    col = tour[1]
    sol[(row-1)*num_nodes+col] = 1
    for i = 1:num_nodes-1
        row = tour[i]
        col = tour[i+1]
        sol[(row-1)*num_nodes+col] = 1
    end
    return sol
end

"""
 Function Name   : Grasp
 Description     : "Greedy ramdomized adaptive search procedure" porpouse is
                    to find a primary solution to a
                    Linear Integer programim problem.

"""
function grasp(m::JuMP.Model,flag_tsp::Bool = true, maxiter::Int = 200, α::Float64 = 0.25, max_iter_improv::Int= 200,max_improvs_aux::Int=200)
    if flag_tsp
        status = _grasp_tsp(m::JuMP.Model,maxiter, α)
    else
        return :NotImplemented
    end
end

"""
 Function Name   : Grasp generic mip problems
 Description     : "Greedy ramdomized adaptive search procedure" porpouse is
                    to find a primary solution to a
                    Linear Integer programim problem.
                    This inplementation was designed to solve any mip problem.
                    (not problem optimized consequently).

"""
function _grasp_generic_mip_problems(model::JuMP.Model,maxiter::Int = 200)
    ## First solution ##
    # Initialization #
    c = []
    variables = collect(1:model.numCols)
    indexVariables = variables[model.colCat .== :Bin]  # list of candidates
    i = 1

    # Fill solution "c" #
    while i < size(indexVariables,1) && i < maxiter
        # Initialization of restricted list of candidates (rlc)
        rlc = []
        for i in size(indexVariables,1)
            # load candidate
            candidate = deepcopy(model)
            candidate.colUpper[indVariable[i]] = 1
            candidate.colLower[indVariable[i]] = 1
            # check feasability
            if _check_feasability_variable(candidate)

            end

        end

        i += 1
    end

end

"""
Function Name   :
Description     :
"""
function _rlc_criteria(model::JuMP.Model)
  indexVariables = collect(1:model.numCols)
  selection = [model.colCat .== :Bin  model.colVal .!= 0  model.colVal .!= 1]
  for row = 1 : size(selection, 1)
    if prod(selection[row, :])
      return indexVariables[row]
    end
  end
  return 0
end

"""
Function Name   : selectInfeasableVariable
Description     : Searches for a varaible that should be binary but has a none
                  binary value in the current model solution.
"""
function _selectInfeasableVariable(model::JuMP.Model)
  indexVariables = collect(1:model.numCols)
  selection = [model.colCat .== :Bin  model.colVal .!= 0  model.colVal .!= 1]
  for row = 1 : size(selection, 1)
    if prod(selection[row, :])
      return indexVariables[row]
    end
  end
  return 0
end

"""
Function Name   : check feasability solution
Description     : Checks if the solution "c" passed as parameter is a feasible
                  solution for the given problem model when linearly relaxed.
"""
function _check_feasability_solution(m::JuMP.Model,c,indVariable)
    for i = 1:size(c,1)
        m.colUpper[indVariable[i]] = c[i]
        m.colLower[indVariable[i]] = c[i]
    end
    obj = m.obj
    obf = JuMP.prepAffObjective(m)
    @objective(m, m.objSense, 0)
    status = solve(m, relaxation=true)
    obval = obf'*m.colVal
    @objective(m, m.objSense, obj)
    if status == :Infeasible
        return false,obval
    end
    return true,obval
end

"""
Function Name   : check feasability variable
Description     : Checks if model has a linearly relaxed feasible solution
"""
function _check_feasability_variable(model::JuMP.Model)
    obj = model.obj
    @objective(m, model.objSense, 0)
    status = solve(model, relaxation=true)
    @objective(m, model.objSense, obj)
    if status == :Infeasible
        return false
    end
    return true
end

"""
Function Name   : Feasability Pump
"""
function feasability_pump()

end

end