text:
module InfiniteLinearAlgebra
using InfiniteArrays: InfRanges
using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, LazyBandedMatrices, SemiseparableMatrices,
        FillArrays, InfiniteArrays, MatrixFactorizations, ArrayLayouts, LinearAlgebra

import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, size, axes, length,
            AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector, Slice,
            show, getproperty, copy, copyto!, map, require_one_based_indexing, similar, inv
import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted

import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout,
                        sublayout, _qr, __qr, MatLmulVec, MatLmulMat, AbstractQLayout, materialize!, diagonaldata, subdiagonaldata, supdiagonaldata,
                        _bidiag_forwardsub!, mulreduce, RangeCumsum, _factorize, transposelayout, ldiv!, lmul!, mul
import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
                        _default_banded_broadcast, banded_similar
import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
import InfiniteArrays: OneToInf, InfUnitRange, Infinity, PosInfinity, InfiniteCardinal, InfStepRange, AbstractInfUnitRange, InfAxes, InfRanges
import LinearAlgebra: matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose, AdjOrTrans
import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, LazyArrayStyle,
                    resizedata!, MemoryLayout,
                    factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
                    applylayout, ApplyLayout, PaddedLayout, zero!, MulAddStyle,
                    LazyArray, LazyMatrix, LazyVector, paddeddata
import MatrixFactorizations: ul, ul!, _ul, ql, ql!, _ql, QLPackedQ, getL, getR, getU, reflector!, reflectorApply!, QL, QR, QRPackedQ,
                            QRPackedQLayout, AdjQRPackedQLayout, QLPackedQLayout, AdjQLPackedQLayout, LayoutQ

import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout, BlockSlice

import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix, banded_chol!

import LazyBandedMatrices: ApplyBandedLayout, BroadcastBandedLayout, _krontrav_axes, _block_interlace_axes, LazyBandedLayout,AbstractLazyBandedLayout

import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMatrix, blockstart, blockstride,
        BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal,
        AbstractBlockBandedLayout, _blockbanded_qr!, BlockBandedLayout

import DSP: conv

import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!

if VERSION < v"1.6-"
    oneto(n) = Base.OneTo(n)
else
    import Base: oneto, unitrange
end

if VERSION ≥ v"1.7-"
    LinearAlgebra._cut_B(x::AbstractVector, r::InfUnitRange) = x
    LinearAlgebra._cut_B(X::AbstractMatrix, r::InfUnitRange) = X
end

# BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()

function ArrayLayouts._power_by_squaring(_, ::NTuple{2,InfiniteCardinal{0}}, A::AbstractMatrix{T}, p::Integer) where T
    if p < 0
        inv(A)^(-p)
    elseif p == 0
        Eye{T}(∞)
    elseif p == 1
        copy(A)
    else
        A*A^(p-1)
    end
end

function choplength(c::AbstractVector, tol)
    @inbounds for k = length(c):-1:1
        if abs(c[k]) > tol
            return k
            break
        end
    end
    return 0
end

_chop!(_, c::AbstractVector, tol::Real) = resize!(c, choplength(c, tol))
_chop!(ax::BlockedUnitRange, c::AbstractVector, tol::Real) = resize!(c, findblock(ax, choplength(c, tol)))

chop!(c::AbstractVector{T}, tol::Real=zero(real(T))) where T = _chop!(axes(c,1), c, tol)

function chop(A::AbstractMatrix{T}, tol::Real=zero(real(T))) where T
    for k = size(A,1):-1:1
        if norm(view(A,k,:))>tol
            A=A[1:k,:]
            break
        end
    end
    for j = size(A,2):-1:1
        if norm(view(A,:,j))>tol
            A=A[:,1:j]
            break
        end
    end
    return A
end

export Vcat, Fill, ql, ql!, ∞, ContinuousSpectrumError, BlockTridiagonal

include("infconv.jl")

include("banded/hessenbergq.jl")

include("banded/infbanded.jl")
include("blockbanded/blockbanded.jl")
include("banded/infqltoeplitz.jl")
include("infql.jl")
include("infqr.jl")
include("inful.jl")
include("infcholesky.jl")

end # module


text:
using KhepriFrame3DD
using Test

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

ext: jl


text:
# File with the function necessary to implement the fast convolution
# in 3D and all the necessary machinery to build the preconditioner

type FastM3D
    ## we may want to add an FFT plan to make the evaluation faster
    # type to encapsulate the fast application of M = I + omega^2G*spadiagm(nu)
    GFFT :: Array{Complex128,3}
    nu :: Array{Float64,1}
    # number of points in the extended domain
    ne :: Int64
    me :: Int64
    le :: Int64
    # number of points in the original domain
    n  :: Int64
    m  :: Int64
    l  :: Int64
    # frequency
    omega :: Float64
    quadRule :: String
    function FastM3D(GFFT,nu,ne,me,le,n,m,l,k;quadRule::String = "Greengard_Vico")
      return new(GFFT,nu,ne,me,le,n,m,l,k,quadRule)
    end
end

import Base.*

function *(M::FastM3D, b::Array{Complex128,1}; verbose::Bool=false)
    # multiply by nu, compute the convolution and then
    # multiply by omega^2
    B = M.omega^2*(FFTconvolution(M,M.nu.*b, verbose=verbose))

    return (b + B)
end

@inline function FFTconvolution(M::FastM3D, b::Array{Complex128,1};
                                verbose::Bool=false )
    # function to compute the convolution with the convolution kernel
    # defined within the FastM3D type using the FFT
    # TODO add a fft plan in here to accelerate the speed
    # input: M::FastM3D type containing the convolution kernel
    #        b::Array{Complex128,1} vector to apply the conv kernel
    verbose && println("Application of the 3D convolution")
    # Allocate the space for the extended B
    BExt = zeros(Complex128,M.ne, M.ne, M.le);
    # zero padding
    BExt[1:M.n,1:M.m,1:M.l]= reshape(b,M.n,M.m,M.l) ;

    # Fourier Transform
    BFft = fftshift(fft(BExt))
    # Component-wise multiplication
    BFft = M.GFFT.*BFft
    # Inverse Fourier Transform
    BExt = ifft(ifftshift(BFft))

    # multiplication by omega^2
    B = BExt[1:M.n,1:M.m,1:M.l];

    return B[:]
end

@everywhere function sampleG3D(k,X,Y,Z, indS, fastconv::FastM3D)
  # function to sample the 3D Green's function in the nodes given by indS
  # input:    k: float64 frequency
  #           X: mesh contaning the X position of each point
  #           Y: mesh contaning the Y position of each point
  #           Z: mesh contaning the Y position of each point
  #           indS: indices in which the sources are located
  #           fastconv: FastM3D type for the application of the
  #                     discrete convolution kernel

  R  = zeros(Complex128, length(indS), length(X))
   for i = 1:length(indS)
      #for i = 1:length(indS)
      ii = indS[i]
      R[i,ii] = 1;
    end

   Gc =  zeros(Complex128, length(indS), length(X))
   # this can be parallelized but then, we may have cache
   # aceess problem
    for i = 1:length(indS)
        Gc[i,:]= FFTconvolution(fastconv, R[i,:][:])
    end
    return Gc
end

@everywhere function sampleG3D(k,X,Y,Z, indS, D0)
  # function to sample the 3D Green's function in the nodes given by indS
  # input:    k: float64 frequency
  #           X: mesh contaning the X position of each point
  #           Y: mesh contaning the Y position of each point
  #           Z: mesh contaning the Y position of each point
  #           indS: indices in which the sources are located
  #           D0: diagonal modification based on the kh, this is given as a
  #               parameter but it should be correctly handled in newer versions
  #               of the code.
  # function to sample the Green's function at frequency k

  R  = SharedArray(Float64, length(indS), length(X))
  Xshared = convert(SharedArray, X)
  Yshared = convert(SharedArray, Y)
  Zshared = convert(SharedArray, Z)
  @sync begin
    @parallel for i = 1:length(indS)
      #for i = 1:length(indS)
      ii = indS[i]
      R[i,:]  = sqrt( (Xshared-Xshared[ii]).^2 + (Yshared-Yshared[ii]).^2 + (Zshared-Zshared[ii]).^2);
      R[i,ii] = 1;
    end
  end

  # sampling the Green's function in the given points
    Gc = (exp(1im*k*R)*h^2)./(4*pi*R) ;
    for i = 1:length(indS)
        ii = indS[i]
        Gc[i,ii]= 1im/4*D0*h^2;
    end
    return Gc
end

#######################################################
### function to enhance parallelism
##################################################
# function SampleGtruncated3D(GFFT, L,k,kx,ky,kz)
#     n = length(kx)
#     m = length(ly)
#     l = length(kz)

#     R  = SharedArray(Float64, length(indS), length(X))
#     kxshared = convert(SharedArray, kx)
#     kyshared = convert(SharedArray, ky)
#     kzshared = convert(SharedArray, kz)

# end

# @everywhere function myrange(q::SharedArray)
#     idx = indexpids(q)
#     if idx == 0
#         # This worker is not assigned a piece
#         return 1:0, 1:0
#     end
#     nchunks = length(procs(q))
#     splits = [round(Int, s) for s in linspace(0,size(q,2),nchunks+1)]
#     1:size(q,1), splits[idx]+1:splits[idx+1]
# end

# @everywhere function sampleGkernelparTruncated3D(k,r::Array{Float64,1},h)
#   n  = length(r)
#   println("Sample kernel parallel loop ")
#   G = SharedArray(Complex128,n)
#   rshared = convert(SharedArray, r)
#   @sync @parallel for ii = 1:n
#           @inbounds  G[ii] = 1im/4*hankelh1(0, k*rshared[ii])*h^2;
#   end
#   return sdata(G)
# end

# ## two different versions of the same function with slight different input

# @everywhere function sampleGkernelpar(k,R::Array{Float64,2},h)
#   (m,n)  = size(R)
#   println("Sample kernel parallel loop with chunks ")
#   G = SharedArray(Complex128,m,n)
#   @time Rshared = convert(SharedArray, R)
#   @sync begin
#         for p in procs(G)
#              @async remotecall_fetch(sampleGkernel_shared_chunk!,p, G, Rshared,k,h)
#         end
#     end
#   return sdata(G)
# end

# # little convenience wrapper
# @everywhere sampleGkernel_shared_chunk!(q,u,k,h) = sampleGkernel_chunk!(q,u,k,h, myrange(q)...)

# @everywhere @inline function sampleGkernel_chunk!(G,R,k::Float64,h::Float64,
#                                                   irange::UnitRange{Int64}, jrange::UnitRange{Int64})
#     #@show (irange, jrange)  # display so we can see what's happening
#     # println(myid())
#     # println(typeof(irange))
#     alpha = 1im/4*h^2
#     for i in irange
#       for j in jrange
#         @inbounds G[i,j] = alpha*hankelh1(0, k*R[i,j]);
#       end
#     end
# end

# we need to write the convolution in 3D, the aim is to have 2 and 3 convolution
function buildFastConvolution3D(x,y,z,X,Y,Z,h,k,nu; quadRule::String = "Greengard_Vico")

 if quadRule == "Greengard_Vico"

      Lp = 4*(abs(x[end] - x[1]) + h)
      L = (abs(x[end] - x[1]) + h)*1.8
      (n,m,l) = length(x), length(y), length(z)
      #LPhysical = abs(x[end]-x[1])+h;

      if mod(n,2) == 0

        kx = (2*pi/Lp)*collect(-(2*n):1:(2*n-1));
        ky = (2*pi/Lp)*collect(-(2*m):1:(2*m-1));
        kz = (2*pi/Lp)*collect(-(2*l):1:(2*l-1));

        # KX = 2*pi*[ kx[i] + 0*j   + 0*p   for i=1:4*n, j=1:4*m, p=1:4*l ]
        # KY = 2*pi*[ 0*i   + ky[j] + 0*p   for i=1:4*n, j=1:4*m, p=1:4*l ]
        # KZ = 2*pi*[ 0*i   + 0*j   + kz[p] for i=1:4*n, j=1:4*m, p=1:4*l ]

        # S = sqrt(KX.^2 + KY.^2 + KZ.^2);
        ## S  = [ sqrt(kx[i]^2 + ky[j]^2 + kz[p]^2) for i=1:4*n, j=1:4*m, p=1:4*l ]

        # GFFT = Gtruncated3D(L, k, S)

        ### GFFT = [ Gtruncated3D(L,k,sqrt(kx[i]^2 + ky[j]^2 + kz[p]^2)) for i=1:4*n, j=1:4*m, p=1:4*l]

        # Computing the convolution kernel (we just use a for loop in order to save memory)
        GFFT = zeros(Complex128, 4*n, 4*m, 4*l)

        for ii=1:4*n, jj=1:4*m, pp=1:4*l
          GFFT[ii,jj,pp]  = Gtruncated3D(L,k,sqrt(kx[ii]^2 + ky[jj]^2 + kz[pp]^2));
        end

        return FastM3D(GFFT,nu(X,Y,Z),4*n,4*m,4*l, n, m,l, k ,quadRule = "Greengard_Vico");
      else

        kx = 2*pi*collect(-2*(n-1):1:2*(n-1) )/4;
        ky = 2*pi*collect(-2*(m-1):1:2*(m-1) )/4;
        kz = 2*pi*collect(-2*(m-1):1:2*(m-1) )/4;

        # KX = [ kx[i] + 0*j   + 0*p   for i=1:4*n-3, j=1:4*m-3, p=1:4*l-3 ]
        # KY = [ 0*i   + ky[j] + 0*p   for i=1:4*n-3, j=1:4*m-3, p=1:4*l-3 ]
        # KZ = [ 0*i   + 0*j   + kz[p] for i=1:4*n-3, j=1:4*m-3, p=1:4*l-3 ]

        # S = sqrt(KX.^2 + KY.^2 + KZ.^2);

        ## S  = [ sqrt(kx[i]^2 + ky[j]^2 + kz[p]^2) for i=1:4*n-3, j=1:4*m-3, p=1:4*l-3 ]

        # GFFT = Gtruncated3D(L, k, S)

        ### GFFT = [ Gtruncated3D(L,k,sqrt(kx[i]^2 + ky[j]^2 + kz[p]^2)) for i=1:4*n-3, j=1:4*m-3, p=1:4*l-3 ]

        # Computing the convolution kernel (we just use a for loop in order to save memory)
        GFFT = zeros(Complex128, 4*n-3, 4*m-3, 4*l-3)

        # This loop can be easily parallelized
        for i=1:4*n-3, j=1:4*m-3, p=1:4*l-3
          GFFT[i,j,p]  = Gtruncated3D(L,k,sqrt(kx[i]^2 + ky[j]^2 + kz[p]^2));
        end

        return FastM3D(GFFT,nu(X,Y,Z),4*n-3,4*m-3,4*l-3,n, m, l,k,quadRule = "Greengard_Vico");

    end
  end
end

function buildSparseAG3DConv(k::Float64,X::Array{Float64,1},Y::Array{Float64,1},
                          Z::Array{Float64,1},
                          fastconv::FastM3D, n::Int64 ,m::Int64,l::Int64; method::String = "normal")
# function that build the sparsigying preconditioner

    # Quick fix: TODO change this!

    Entries = entriesSparseG3D(k,X,Y,Z,fastconv, n ,m, l);
    Ind = reshape(collect(1:n*m*l),n,m,l);

    if method == "normal"
      (Indices, Values) = entriesSparseA3D(k,X,Y,Z,fastconv, n ,m,l);
    elseif method == "randomized"
      println("Not implemented yet")
      (Indices, Values) = entriesSparseARand(k,X,Y,D0, n ,m);
    end

    ValuesAG = Values[1]*Entries[1];

    # building the indices, columns and rows for the interior
    (rowA, colA, valA) = createIndices(Ind[2:end-1,2:end-1,2:end-1][:],
                                    Indices[1][:], (Values[1]*Entries[1])[ :]);

    # for  x = xmin,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,2:end-1][:],
                                    Indices[2][:], (Values[2]*Entries[2])[ :]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmax,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,2:end-1][:],
                                    Indices[3][:], (Values[3]*Entries[3])[ :]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  y = ymin,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,2:end-1][:],
                                    Indices[4][:], (Values[4]*Entries[4])[ :]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  y = ymax,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,2:end-1][:],
                                    Indices[5][:], (Values[5]*Entries[5])[ :]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  z = zmin,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,1][:],
                                    Indices[6][:], (Values[6]*Entries[6])[ :]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  z = zmax,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,end][:],
                                    Indices[7][:], (Values[7]*Entries[7])[ :]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

# we need to incorporate the vertices

# indvertex1  = round(Integer, changeInd3D(1,1,lHalf,n,m,l) + subStencil3D(3,3,2,Ind_relative));
# for  x = xmin,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[1,1,2:end-1][:],
                                    Indices[8][:], (Values[8]*Entries[8])[ :]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex2  = round(Integer, changeInd3D(n,1,lHalf,n,m,l) + subStencil3D(1,3,2,Ind_relative));
# for  x = xmax,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[end,1,2:end-1][:],
                                    Indices[9][:], (Values[9]*Entries[9])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex3  = round(Integer, changeInd3D(1,m,lHalf,n,m,l) + subStencil3D(3,1,2,Ind_relative));
# for  x = xmin,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,end,2:end-1][:],
                                    Indices[10][:], (Values[10]*Entries[10])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex4  = round(Integer, changeInd3D(n,m,lHalf,n,m,l) + subStencil3D(1,1,2,Ind_relative));
# for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[end,end,2:end-1][:],
                                    Indices[11][:], (Values[11]*Entries[11])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex5  = round(Integer, changeInd3D(1,mHalf,1,n,m,l) + subStencil3D(3,2,3,Ind_relative));
# for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,1][:],
                                    Indices[12][:], (Values[12]*Entries[12])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex6  = round(Integer, changeInd3D(n,mHalf,1,n,m,l) + subStencil3D(1,2,3,Ind_relative));
# for  x = xmax,  y = any z = zmin
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,1][:],
                                    Indices[13][:], (Values[13]*Entries[13])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex7  = round(Integer, changeInd3D(1,mHalf,l,n,m,l) + subStencil3D(3,2,1,Ind_relative));
# for  x = xmin,  y = any z = zmax
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,end][:],
                                    Indices[14][:], (Values[14]*Entries[14])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex8  = round(Integer, changeInd3D(n,mHalf,l,n,m,l) + subStencil3D(1,2,1,Ind_relative));
# for  x = xmax,  y = any  z = zmax
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,end][:],
                                    Indices[15][:], (Values[15]*Entries[15])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex9  = round(Integer, changeInd3D(nHalf,1,1,n,m,l) + subStencil3D(2,3,3,Ind_relative));
# for  x = any,  y = tmin  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,1][:],
                                    Indices[16][:], (Values[16]*Entries[16])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex10 = round(Integer, changeInd3D(nHalf,m,1,n,m,l) + subStencil3D(2,1,3,Ind_relative));
# for  x = any,  y = ymax  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,1][:],
                                    Indices[17][:], (Values[17]*Entries[17])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex11 = round(Integer, changeInd3D(nHalf,1,l,n,m,l) + subStencil3D(2,3,1,Ind_relative));
# for  x = any,  y = ymin  z = zmax
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,end][:],
                                    Indices[18][:], (Values[18]*Entries[18])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex12 = round(Integer, changeInd3D(nHalf,m,l,n,m,l) + subStencil3D(2,1,1,Ind_relative));
# for  x = any,  y = ymax  z = zmax
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,end][:],
                                    Indices[19][:], (Values[19]*Entries[19])[:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

#################################################################################
############## Now we incorporate the corners  ##################################
#################################################################################

# indcorner1 = round(Integer, changeInd3D(1,1,1,n,m,l) + Ind_relative[2:3,2:3,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[1,1,1],
                                    Indices[20][:], (Values[20]*Entries[20])[:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner2 = round(Integer, changeInd3D(n,1,1,n,m,l) + Ind_relative[1:2,2:3,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[end,1,1],
                                    Indices[21][:], (Values[21]*Entries[21])[:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner3 = round(Integer, changeInd3D(1,m,1,n,m,l) + Ind_relative[2:3,1:2,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[1,end,1],
                                    Indices[22][:], (Values[22]*Entries[22])[:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner4 = round(Integer, changeInd3D(n,m,1,n,m,l) + Ind_relative[1:2,1:2,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[end,end,1],
                                    Indices[23][:], (Values[23]*Entries[23])[:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner5 = round(Integer, changeInd3D(1,1,l,n,m,l) + Ind_relative[2:3,2:3,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[1,1,end],
                                    Indices[24][:], (Values[24]*Entries[24])[:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner6 = round(Integer, changeInd3D(n,1,l,n,m,l) + Ind_relative[1:2,2:3,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[end,1,end],
                                    Indices[25][:], (Values[25]*Entries[25])[:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner7 = round(Integer, changeInd3D(1,m,l,n,m,l) + Ind_relative[2:3,1:2,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[1,end,end],
                                    Indices[26][:], (Values[26]*Entries[26])[:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner8 = round(Integer, changeInd3D(n,m,l,n,m,l) + Ind_relative[1:2,1:2,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[end,end,end],
                                    Indices[27][:], (Values[27]*Entries[27])[:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    AG = sparse(rowA,colA,valA);

    return AG;
end

@everywhere function sampleG3D(k,X,Y,Z, indS, fastconv::FastM3D)
  # function to sample the 3D Green's function in the nodes given by indS
  # input:    k: float64 frequency
  #           X: mesh contaning the X position of each point
  #           Y: mesh contaning the Y position of each point
  #           Z: mesh contaning the Y position of each point
  #           indS: indices in which the sources are located
  #           D0: how to compute the qudrature fast
  # function to sample the Green's function at frequency k

  R  = zeros(Complex128, length(indS), length(X))
   for i = 1:length(indS)
      #for i = 1:length(indS)
      ii = indS[i]
      R[i,ii] = 1;
    end

   Gc =  zeros(Complex128, length(indS), length(X))
    for i = 1:length(indS)
        Gc[i,:]= FFTconvolution(fastconv, R[i,:][:])
    end
    return Gc
end

function entriesSparseG3D(k::Float64,X::Array{Float64,1},Y::Array{Float64,1},
                          Z::Array{Float64,1},
                          fastconv::FastM3D, n::Int64 ,m::Int64,l::Int64)
  # function to compute the entried of G, inside the volume, at the boundaries
  # and at the corners. This allows us to compute A*G in O(n) time instead of
  # O(n^2)
  # we need to have an even number of points
  #

  Entries  = Array{Complex128}[]

  N = n*m*l;

  Ind_relative = zeros(Int64,3,3,3)
  Ind_relative[:,:,1] = [(-m*n-n-1) (-m*n-n) (-m*n-n+1);
                         (-m*n  -1) (-m*n  ) (-m*n  +1);
                         (-m*n+n-1) (-m*n+n) (-m*n+n+1) ]';

  Ind_relative[:,:,2] = [-n-1 -n -n+1;
                           -1  0    1;
                          n-1  n  n+1]';

  Ind_relative[:,:,3] = [(m*n-n-1) (m*n-n) (m*n-n+1);
                         (m*n  -1) (m*n  ) (m*n  +1);
                         (m*n+n-1) (m*n+n) (m*n+n+1) ]';

  # computing the entries for the interior

  nHalf = round(Integer,n/2);
  mHalf = round(Integer,m/2);
  lHalf = round(Integer,l/2);

  indVol = round(Integer, changeInd3D(nHalf,mHalf,lHalf,n,m,l)+Ind_relative[:]);

  GSampled = sampleG3D(k,X,Y,Z,indVol, fastconv)[:,indVol];

  push!(Entries,GSampled);

  # for  x = xmin,  y = anything z = anything
  indFx1 = round(Integer, changeInd3D(1,mHalf,lHalf,n,m,l) + Ind_relative[2:3,:,:][:] );
  GSampled = sampleG3D(k,X,Y,Z,indFx1, fastconv)[:,indFx1];

  push!(Entries,GSampled);

  # for  x = xmax, y = any z = any
  indFxN = round(Integer, changeInd3D(n,mHalf,lHalf,n,m,l) + Ind_relative[1:2,:,:][:]);
  GSampled = sampleG3D(k,X,Y,Z,indFxN, fastconv)[:,indFxN];

  push!(Entries,GSampled);

  # for  y = ymin, x = any z = any
  indFy1 = round(Integer, changeInd3D(nHalf,1,lHalf,n,m,l) + Ind_relative[:,2:3,:][:] );
  GSampled = sampleG3D(k,X,Y,Z,indFy1, fastconv)[:,indFy1];

  push!(Entries,GSampled);

  # for  y = ymax, x = any z = any
  indFyN = round(Integer, changeInd3D(nHalf,m,lHalf,n,m,l) + Ind_relative[:,1:2,:][:] );
  GSampled = sampleG3D(k,X,Y,Z,indFyN, fastconv)[:,indFyN];

  push!(Entries,GSampled);

  # for  z = zmin, x = any y = any
  indFz1 = round(Integer, changeInd3D(nHalf,mHalf,1,n,m,l) + Ind_relative[:,:,2:3][:] );
  GSampled = sampleG3D(k,X,Y,Z,indFz1, fastconv)[:,indFz1];

  push!(Entries,GSampled);

  # for  z = zmax, x = any y = any
  indFzN = round(Integer, changeInd3D(nHalf,mHalf,l,n,m,l) + Ind_relative[:,:,1:2][:] );
  GSampled = sampleG3D(k,X,Y,Z,indFzN, fastconv)[:,indFzN];

  push!(Entries,GSampled);

  # we need to incorporate the vertices
  indvertex1  = round(Integer, changeInd3D(1,1,lHalf,n,m,l) + subStencil3D(3,3,2,Ind_relative));
  indvertex2  = round(Integer, changeInd3D(n,1,lHalf,n,m,l) + subStencil3D(1,3,2,Ind_relative));
  indvertex3  = round(Integer, changeInd3D(1,m,lHalf,n,m,l) + subStencil3D(3,1,2,Ind_relative));
  indvertex4  = round(Integer, changeInd3D(n,m,lHalf,n,m,l) + subStencil3D(1,1,2,Ind_relative));
  indvertex5  = round(Integer, changeInd3D(1,mHalf,1,n,m,l) + subStencil3D(3,2,3,Ind_relative));
  indvertex6  = round(Integer, changeInd3D(n,mHalf,1,n,m,l) + subStencil3D(1,2,3,Ind_relative));
  indvertex7  = round(Integer, changeInd3D(1,mHalf,l,n,m,l) + subStencil3D(3,2,1,Ind_relative));
  indvertex8  = round(Integer, changeInd3D(n,mHalf,l,n,m,l) + subStencil3D(1,2,1,Ind_relative));
  indvertex9  = round(Integer, changeInd3D(nHalf,1,1,n,m,l) + subStencil3D(2,3,3,Ind_relative));
  indvertex10 = round(Integer, changeInd3D(nHalf,m,1,n,m,l) + subStencil3D(2,1,3,Ind_relative));
  indvertex11 = round(Integer, changeInd3D(nHalf,1,l,n,m,l) + subStencil3D(2,3,1,Ind_relative));
  indvertex12 = round(Integer, changeInd3D(nHalf,m,l,n,m,l) + subStencil3D(2,1,1,Ind_relative));

  GSampled = sampleG3D(k,X,Y,Z,indvertex1, fastconv)[:,indvertex1 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex2, fastconv)[:,indvertex2 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex3, fastconv)[:,indvertex3 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex4, fastconv)[:,indvertex4 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex5, fastconv)[:,indvertex5 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex6, fastconv)[:,indvertex6 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex7, fastconv)[:,indvertex7 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex8, fastconv)[:,indvertex8 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex9, fastconv)[:,indvertex9 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex10, fastconv)[:,indvertex10 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex11, fastconv)[:,indvertex11 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indvertex12, fastconv)[:,indvertex12 ];
  push!(Entries,GSampled);

  # Now we incorporate the corners
  indcorner1 = round(Integer, changeInd3D(1,1,1,n,m,l) + Ind_relative[2:3,2:3,2:3][:]);
  indcorner2 = round(Integer, changeInd3D(n,1,1,n,m,l) + Ind_relative[1:2,2:3,2:3][:]);
  indcorner3 = round(Integer, changeInd3D(1,m,1,n,m,l) + Ind_relative[2:3,1:2,2:3][:]);
  indcorner4 = round(Integer, changeInd3D(n,m,1,n,m,l) + Ind_relative[1:2,1:2,2:3][:]);
  indcorner5 = round(Integer, changeInd3D(1,1,l,n,m,l) + Ind_relative[2:3,2:3,1:2][:]);
  indcorner6 = round(Integer, changeInd3D(n,1,l,n,m,l) + Ind_relative[1:2,2:3,1:2][:]);
  indcorner7 = round(Integer, changeInd3D(1,m,l,n,m,l) + Ind_relative[2:3,1:2,1:2][:]);
  indcorner8 = round(Integer, changeInd3D(n,m,l,n,m,l) + Ind_relative[1:2,1:2,1:2][:]);

  GSampled = sampleG3D(k,X,Y,Z,indcorner1, fastconv)[:,indcorner1 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indcorner2, fastconv)[:,indcorner2 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indcorner3, fastconv)[:,indcorner3 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indcorner4, fastconv)[:,indcorner4 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indcorner5, fastconv)[:,indcorner5 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indcorner6, fastconv)[:,indcorner6 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indcorner7, fastconv)[:,indcorner7 ];
  push!(Entries,GSampled);

  GSampled = sampleG3D(k,X,Y,Z,indcorner8, fastconv)[:,indcorner8 ];
  push!(Entries,GSampled);

  return Entries

end

function entriesSparseA3D(k,X::Array{Float64,1},Y::Array{Float64,1},
                          Z::Array{Float64,1},fastconv::FastM3D,
                          n::Int64 ,m::Int64, l::Int64)
  # in this case we need to build everythig with ranodmized methods
  # we need to have an odd number of points
  #@assert mod(length(X),2) == 1
  Entries  = Array{Complex128}[]
  Indices  = Array{Int64}[]

  N = n*m*l;

  Ind_relative = zeros(Int64,3,3,3)
  Ind_relative[:,:,1] = [(-m*n-n-1) (-m*n-n) (-m*n-n+1);
                         (-m*n  -1) (-m*n  ) (-m*n  +1);
                         (-m*n+n-1) (-m*n+n) (-m*n+n+1) ]';

  Ind_relative[:,:,2] = [-n-1 -n -n+1;
                           -1  0    1;
                          n-1  n  n+1]';

  Ind_relative[:,:,3] = [(m*n-n-1) (m*n-n) (m*n-n+1);
                         (m*n  -1) (m*n  ) (m*n  +1);
                         (m*n+n-1) (m*n+n) (m*n+n+1) ]';

  # computing the entries for the interior

  nHalf = round(Integer,n/2);
  mHalf = round(Integer,m/2);
  lHalf = round(Integer,l/2);

  indVol = round(Integer, changeInd3D(nHalf,mHalf,lHalf,n,m,l)+Ind_relative[:]);

  indVolC = setdiff(collect(1:N),indVol);
  GSampled = sampleG3D(k,X,Y,Z,indVol, fastconv)[:,indVolC ];

  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative);

  # for  x = xmin,  y = anything z = anything
  indFx1 = round(Integer, changeInd3D(1,mHalf,lHalf,n,m,l) + Ind_relative[2:3,:,:][:] );
  indC = setdiff(collect(1:N),indFx1);
  GSampled = sampleG3D(k,X,Y,Z,indFx1, fastconv)[:,indC ];
   # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  # saving the entries
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[2:3,:,:][:] );#'

  # for  x = xmax, y = any z = any
  indFxN = round(Integer, changeInd3D(n,mHalf,lHalf,n,m,l) + Ind_relative[1:2,:,:][:]);
  indC = setdiff(collect(1:N),indFxN);
  GSampled = sampleG3D(k,X,Y,Z,indFxN, fastconv)[:,indC ];
  (U,s,V) = svd(GSampled);
  # saving the entries
  push!(Entries,U[:,end]'); #'
  push!(Indices, Ind_relative[1:2,:,:][:]); #'

  # for  y = ymin, x = any z = any
  indFy1 = round(Integer, changeInd3D(nHalf,1,lHalf,n,m,l) + Ind_relative[:,2:3,:][:] );
  indC = setdiff(collect(1:N),indFy1);
  GSampled = sampleG3D(k,X,Y,Z,indFy1, fastconv)[:,indC ];
  (U,s,V) = svd(GSampled);
  # saving the entries
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[:,2:3,:][:]); #'

  # for  y = ymax, x = any z = any
  indFyN = round(Integer, changeInd3D(nHalf,m,lHalf,n,m,l) + Ind_relative[:,1:2,:][:] );
  indC = setdiff(collect(1:N),indFyN);
  GSampled = sampleG3D(k,X,Y,Z,indFyN, fastconv)[:,indC ];
  (U,s,V) = svd(GSampled);
  # saving the entries
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[:,1:2,:][:] ); #'

  # for  z = zmin, x = any y = any
  indFz1 = round(Integer, changeInd3D(nHalf,mHalf,1,n,m,l) + Ind_relative[:,:,2:3][:] );
  indC = setdiff(collect(1:N),indFz1);
  GSampled = sampleG3D(k,X,Y,Z,indFz1, fastconv)[:,indC ];
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[:,:,2:3][:] ); #'

  # for  z = zmax, x = any y = any
  indFzN = round(Integer, changeInd3D(nHalf,mHalf,l,n,m,l) + Ind_relative[:,:,1:2][:] );
  indC = setdiff(collect(1:N),indFzN);
  GSampled = sampleG3D(k,X,Y,Z,indFzN, fastconv)[:,indC ];
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[:,:,1:2][:] ); #'

  # we need to incorporate the vertices
  indvertex1  = round(Integer, changeInd3D(1,1,lHalf,n,m,l) + subStencil3D(3,3,2,Ind_relative));
  indvertex2  = round(Integer, changeInd3D(n,1,lHalf,n,m,l) + subStencil3D(1,3,2,Ind_relative));
  indvertex3  = round(Integer, changeInd3D(1,m,lHalf,n,m,l) + subStencil3D(3,1,2,Ind_relative));
  indvertex4  = round(Integer, changeInd3D(n,m,lHalf,n,m,l) + subStencil3D(1,1,2,Ind_relative));
  indvertex5  = round(Integer, changeInd3D(1,mHalf,1,n,m,l) + subStencil3D(3,2,3,Ind_relative));
  indvertex6  = round(Integer, changeInd3D(n,mHalf,1,n,m,l) + subStencil3D(1,2,3,Ind_relative));
  indvertex7  = round(Integer, changeInd3D(1,mHalf,l,n,m,l) + subStencil3D(3,2,1,Ind_relative));
  indvertex8  = round(Integer, changeInd3D(n,mHalf,l,n,m,l) + subStencil3D(1,2,1,Ind_relative));
  indvertex9  = round(Integer, changeInd3D(nHalf,1,1,n,m,l) + subStencil3D(2,3,3,Ind_relative));
  indvertex10 = round(Integer, changeInd3D(nHalf,m,1,n,m,l) + subStencil3D(2,1,3,Ind_relative));
  indvertex11 = round(Integer, changeInd3D(nHalf,1,l,n,m,l) + subStencil3D(2,3,1,Ind_relative));
  indvertex12 = round(Integer, changeInd3D(nHalf,m,l,n,m,l) + subStencil3D(2,1,1,Ind_relative));

  indC = setdiff(collect(1:N),indvertex1);
  GSampled = sampleG3D(k,X,Y,Z,indvertex1, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(3,3,2,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex2);
  GSampled = sampleG3D(k,X,Y,Z,indvertex2, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(1,3,2,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex3);
  GSampled = sampleG3D(k,X,Y,Z,indvertex3, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(3,1,2,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex4);
  GSampled = sampleG3D(k,X,Y,Z,indvertex4, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(1,1,2,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex5);
  GSampled = sampleG3D(k,X,Y,Z,indvertex5, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(3,2,3,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex6);
  GSampled = sampleG3D(k,X,Y,Z,indvertex6, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(1,2,3,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex7);
  GSampled = sampleG3D(k,X,Y,Z,indvertex7, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(3,2,1,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex8);
  GSampled = sampleG3D(k,X,Y,Z,indvertex8, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(1,2,1,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex9);
  GSampled = sampleG3D(k,X,Y,Z,indvertex9, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(2,3,3,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex10);
  GSampled = sampleG3D(k,X,Y,Z,indvertex10, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(2,1,3,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex11);
  GSampled = sampleG3D(k,X,Y,Z,indvertex11, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(2,3,1,Ind_relative)); #'

  indC = setdiff(collect(1:N),indvertex12);
  GSampled = sampleG3D(k,X,Y,Z,indvertex12, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,subStencil3D(2,1,1,Ind_relative)); #'

  # Now we incorporate the corners
  indcorner1 = round(Integer, changeInd3D(1,1,1,n,m,l) + Ind_relative[2:3,2:3,2:3][:]);
  indcorner2 = round(Integer, changeInd3D(n,1,1,n,m,l) + Ind_relative[1:2,2:3,2:3][:]);
  indcorner3 = round(Integer, changeInd3D(1,m,1,n,m,l) + Ind_relative[2:3,1:2,2:3][:]);
  indcorner4 = round(Integer, changeInd3D(n,m,1,n,m,l) + Ind_relative[1:2,1:2,2:3][:]);
  indcorner5 = round(Integer, changeInd3D(1,1,l,n,m,l) + Ind_relative[2:3,2:3,1:2][:]);
  indcorner6 = round(Integer, changeInd3D(n,1,l,n,m,l) + Ind_relative[1:2,2:3,1:2][:]);
  indcorner7 = round(Integer, changeInd3D(1,m,l,n,m,l) + Ind_relative[2:3,1:2,1:2][:]);
  indcorner8 = round(Integer, changeInd3D(n,m,l,n,m,l) + Ind_relative[1:2,1:2,1:2][:]);

  indC = setdiff(collect(1:N),indcorner1);
  GSampled = sampleG3D(k,X,Y,Z,indcorner1, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[2:3,2:3,2:3][:]); #'

  #'
  indC = setdiff(collect(1:N),indcorner2);
  GSampled = sampleG3D(k,X,Y,Z,indcorner2, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[1:2,2:3,2:3][:]); #'

  #'
  indC = setdiff(collect(1:N),indcorner3);
  GSampled = sampleG3D(k,X,Y,Z,indcorner3, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[2:3,1:2,2:3][:]); #'

  #'
  indC = setdiff(collect(1:N),indcorner4);
  GSampled = sampleG3D(k,X,Y,Z,indcorner4, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[1:2,1:2,2:3][:]); #'

  indC = setdiff(collect(1:N),indcorner5);
  GSampled = sampleG3D(k,X,Y,Z,indcorner5, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[2:3,2:3,1:2][:]); #'

  #'
  indC = setdiff(collect(1:N),indcorner6);
  GSampled = sampleG3D(k,X,Y,Z,indcorner6, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[1:2,2:3,1:2][:]); #'

  #'
  indC = setdiff(collect(1:N),indcorner7);
  GSampled = sampleG3D(k,X,Y,Z,indcorner7, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[2:3,1:2,1:2][:]); #'

  #'
  indC = setdiff(collect(1:N),indcorner8);
  GSampled = sampleG3D(k,X,Y,Z,indcorner8, fastconv)[:,indC ];
  # computing the sparsifying correction
  (U,s,V) = svd(GSampled);
  push!(Entries,U[:,end]'); #'
  push!(Indices,Ind_relative[1:2,1:2,1:2][:]); #'

  return (Indices, Entries)
end

## To be done! (I don't remember how I built this one)
function buildSparseA3DConv(k::Float64,X::Array{Float64,1},Y::Array{Float64,1},
                        Z::Array{Float64,1},
                        fastconv::FastM3D, n::Int64 ,m::Int64,l::Int64; method::String = "normal")
# function that build the sparsigying preconditioner
# this varians uses the kernel from the discretization of the integral system

    Ind = reshape(collect(1:n*m*l),n,m,l);

    if method == "normal"
      (Indices, Values) = entriesSparseA3D(k,X,Y,Z,fastconv, n ,m,l);
    elseif method == "randomized"
      println("Not implemented yet")
      (Indices, Values) = entriesSparseARand(k,X,Y,fastconv, n ,m);
    end

    # building the indices, columns and rows for the interior
    (rowA, colA, valA) = createIndices(Ind[2:end-1,2:end-1,2:end-1][:],
                                    Indices[1][:], Values[1][:]);

    # for  x = xmin,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,2:end-1][:],
                                    Indices[2][:], Values[2][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmax,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,2:end-1][:],
                                    Indices[3][:], Values[3][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  y = ymin,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,2:end-1][:],
                                    Indices[4][:], Values[4][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  y = ymax,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,2:end-1][:],
                                    Indices[5][:], Values[5][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  z = zmin,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,1][:],
                                    Indices[6][:], Values[6][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  z = zmax,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,end][:],
                                    Indices[7][:], Values[7][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

# we need to incorporate the vertices

# indvertex1  = round(Integer, changeInd3D(1,1,lHalf,n,m,l) + subStencil3D(3,3,2,Ind_relative));
# for  x = xmin,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[1,1,2:end-1][:],
                                    Indices[8][:], Values[8][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex2  = round(Integer, changeInd3D(n,1,lHalf,n,m,l) + subStencil3D(1,3,2,Ind_relative));
# for  x = xmax,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[end,1,2:end-1][:],
                                    Indices[9][:], Values[9][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex3  = round(Integer, changeInd3D(1,m,lHalf,n,m,l) + subStencil3D(3,1,2,Ind_relative));
# for  x = xmin,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,end,2:end-1][:],
                                    Indices[10][:], Values[10][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex4  = round(Integer, changeInd3D(n,m,lHalf,n,m,l) + subStencil3D(1,1,2,Ind_relative));
# for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[end,end,2:end-1][:],
                                    Indices[11][:], Values[11][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex5  = round(Integer, changeInd3D(1,mHalf,1,n,m,l) + subStencil3D(3,2,3,Ind_relative));
# for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,1][:],
                                    Indices[12][:], Values[12][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex6  = round(Integer, changeInd3D(n,mHalf,1,n,m,l) + subStencil3D(1,2,3,Ind_relative));
# for  x = xmax,  y = any z = zmin
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,1][:],
                                    Indices[13][:], Values[13][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex7  = round(Integer, changeInd3D(1,mHalf,l,n,m,l) + subStencil3D(3,2,1,Ind_relative));
# for  x = xmin,  y = any z = zmax
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,end][:],
                                    Indices[14][:], Values[14][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex8  = round(Integer, changeInd3D(n,mHalf,l,n,m,l) + subStencil3D(1,2,1,Ind_relative));
# for  x = xmax,  y = any  z = zmax
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,end][:],
                                    Indices[15][:], Values[15][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex9  = round(Integer, changeInd3D(nHalf,1,1,n,m,l) + subStencil3D(2,3,3,Ind_relative));
# for  x = any,  y = tmin  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,1][:],
                                    Indices[16][:], Values[16][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex10 = round(Integer, changeInd3D(nHalf,m,1,n,m,l) + subStencil3D(2,1,3,Ind_relative));
# for  x = any,  y = ymax  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,1][:],
                                    Indices[17][:], Values[17][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex11 = round(Integer, changeInd3D(nHalf,1,l,n,m,l) + subStencil3D(2,3,1,Ind_relative));
# for  x = any,  y = ymin  z = zmax
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,end][:],
                                    Indices[18][:], Values[18][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

# indvertex12 = round(Integer, changeInd3D(nHalf,m,l,n,m,l) + subStencil3D(2,1,1,Ind_relative));
# for  x = any,  y = ymax  z = zmax
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,end][:],
                                    Indices[19][:], Values[19][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

#################################################################################
############## Now we incorporate the corners  ##################################
#################################################################################

# indcorner1 = round(Integer, changeInd3D(1,1,1,n,m,l) + Ind_relative[2:3,2:3,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[1,1,1],
                                    Indices[20][:], Values[20][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner2 = round(Integer, changeInd3D(n,1,1,n,m,l) + Ind_relative[1:2,2:3,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[end,1,1],
                                    Indices[21][:], Values[21][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner3 = round(Integer, changeInd3D(1,m,1,n,m,l) + Ind_relative[2:3,1:2,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[1,end,1],
                                    Indices[22][:], Values[22][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner4 = round(Integer, changeInd3D(n,m,1,n,m,l) + Ind_relative[1:2,1:2,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[end,end,1],
                                    Indices[23][:], Values[23][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner5 = round(Integer, changeInd3D(1,1,l,n,m,l) + Ind_relative[2:3,2:3,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[1,1,end],
                                    Indices[24][:], Values[24][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner6 = round(Integer, changeInd3D(n,1,l,n,m,l) + Ind_relative[1:2,2:3,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[end,1,end],
                                    Indices[25][:], Values[25][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner7 = round(Integer, changeInd3D(1,m,l,n,m,l) + Ind_relative[2:3,1:2,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[1,end,end],
                                    Indices[26][:], Values[26][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)
# indcorner8 = round(Integer, changeInd3D(n,m,l,n,m,l) + Ind_relative[1:2,1:2,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[end,end,end],
                                    Indices[27][:], Values[27][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    A = sparse(rowA,colA,valA);

    return A;
end

function buildSparseAG3DConv(k::Float64,X::Array{Float64,1},Y::Array{Float64,1},
                          Z::Array{Float64,1},
                          fastconv::FastM3D, n::Int64 ,m::Int64,l::Int64; method::String = "normal")
# function that build the sparsigying preconditioner

    # Quick fix: TODO change this!

    Entries = entriesSparseG3D(k,X,Y,Z,fastconv, n ,m, l);
    Ind = reshape(collect(1:n*m*l),n,m,l);

    if method == "normal"
      (Indices, Values) = entriesSparseA3D(k,X,Y,Z,fastconv, n ,m,l);
    elseif method == "randomized"
      println("Not implemented yet")
      (Indices, Values) = entriesSparseARand(k,X,Y,D0, n ,m);
    end

    ValuesAG = Values[1]*Entries[1];

    # building the indices, columns and rows for the interior
    (rowAG, colAG, valAG) = createIndices(Ind[2:end-1,2:end-1,2:end-1][:],
                                    Indices[1][:], (Values[1]*Entries[1])[ :]);

    # for  x = xmin,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,2:end-1][:],
                                    Indices[2][:], (Values[2]*Entries[2])[ :]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

    # for  x = xmax,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,2:end-1][:],
                                    Indices[3][:], (Values[3]*Entries[3])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    # for  y = ymin,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,2:end-1][:],
                                    Indices[4][:], (Values[4]*Entries[4])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    # for  y = ymax,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,2:end-1][:],
                                    Indices[5][:], (Values[5]*Entries[5])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    # for  z = zmin,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,1][:],
                                    Indices[6][:], (Values[6]*Entries[6])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    # for  z = zmax,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,end][:],
                                    Indices[7][:], (Values[7]*Entries[7])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

# we need to incorporate the vertices

# indvertex1  = round(Integer, changeInd3D(1,1,lHalf,n,m,l) + subStencil3D(3,3,2,Ind_relative));
# for  x = xmin,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[1,1,2:end-1][:],
                                    Indices[8][:], (Values[8]*Entries[8])[ :]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex2  = round(Integer, changeInd3D(n,1,lHalf,n,m,l) + subStencil3D(1,3,2,Ind_relative));
# for  x = xmax,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[end,1,2:end-1][:],
                                    Indices[9][:], (Values[9]*Entries[9])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex3  = round(Integer, changeInd3D(1,m,lHalf,n,m,l) + subStencil3D(3,1,2,Ind_relative));
# for  x = xmin,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,end,2:end-1][:],
                                    Indices[10][:], (Values[10]*Entries[10])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex4  = round(Integer, changeInd3D(n,m,lHalf,n,m,l) + subStencil3D(1,1,2,Ind_relative));
# for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[end,end,2:end-1][:],
                                    Indices[11][:], (Values[11]*Entries[11])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex5  = round(Integer, changeInd3D(1,mHalf,1,n,m,l) + subStencil3D(3,2,3,Ind_relative));
# for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,1][:],
                                    Indices[12][:], (Values[12]*Entries[12])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex6  = round(Integer, changeInd3D(n,mHalf,1,n,m,l) + subStencil3D(1,2,3,Ind_relative));
# for  x = xmax,  y = any z = zmin
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,1][:],
                                    Indices[13][:], (Values[13]*Entries[13])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex7  = round(Integer, changeInd3D(1,mHalf,l,n,m,l) + subStencil3D(3,2,1,Ind_relative));
# for  x = xmin,  y = any z = zmax
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,end][:],
                                    Indices[14][:], (Values[14]*Entries[14])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex8  = round(Integer, changeInd3D(n,mHalf,l,n,m,l) + subStencil3D(1,2,1,Ind_relative));
# for  x = xmax,  y = any  z = zmax
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,end][:],
                                    Indices[15][:], (Values[15]*Entries[15])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex9  = round(Integer, changeInd3D(nHalf,1,1,n,m,l) + subStencil3D(2,3,3,Ind_relative));
# for  x = any,  y = tmin  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,1][:],
                                    Indices[16][:], (Values[16]*Entries[16])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex10 = round(Integer, changeInd3D(nHalf,m,1,n,m,l) + subStencil3D(2,1,3,Ind_relative));
# for  x = any,  y = ymax  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,1][:],
                                    Indices[17][:], (Values[17]*Entries[17])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex11 = round(Integer, changeInd3D(nHalf,1,l,n,m,l) + subStencil3D(2,3,1,Ind_relative));
# for  x = any,  y = ymin  z = zmax
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,end][:],
                                    Indices[18][:], (Values[18]*Entries[18])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex12 = round(Integer, changeInd3D(nHalf,m,l,n,m,l) + subStencil3D(2,1,1,Ind_relative));
# for  x = any,  y = ymax  z = zmax
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,end][:],
                                    Indices[19][:], (Values[19]*Entries[19])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

#################################################################################
############## Now we incorporate the corners  ##################################
#################################################################################

# indcorner1 = round(Integer, changeInd3D(1,1,1,n,m,l) + Ind_relative[2:3,2:3,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[1,1,1],
                                    Indices[20][:], (Values[20]*Entries[20])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner2 = round(Integer, changeInd3D(n,1,1,n,m,l) + Ind_relative[1:2,2:3,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[end,1,1],
                                    Indices[21][:], (Values[21]*Entries[21])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner3 = round(Integer, changeInd3D(1,m,1,n,m,l) + Ind_relative[2:3,1:2,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[1,end,1],
                                    Indices[22][:], (Values[22]*Entries[22])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner4 = round(Integer, changeInd3D(n,m,1,n,m,l) + Ind_relative[1:2,1:2,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[end,end,1],
                                    Indices[23][:], (Values[23]*Entries[23])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner5 = round(Integer, changeInd3D(1,1,l,n,m,l) + Ind_relative[2:3,2:3,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[1,1,end],
                                    Indices[24][:], (Values[24]*Entries[24])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner6 = round(Integer, changeInd3D(n,1,l,n,m,l) + Ind_relative[1:2,2:3,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[end,1,end],
                                    Indices[25][:], (Values[25]*Entries[25])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner7 = round(Integer, changeInd3D(1,m,l,n,m,l) + Ind_relative[2:3,1:2,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[1,end,end],
                                    Indices[26][:], (Values[26]*Entries[26])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner8 = round(Integer, changeInd3D(n,m,l,n,m,l) + Ind_relative[1:2,1:2,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[end,end,end],
                                    Indices[27][:], (Values[27]*Entries[27])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    AG = sparse(rowAG,colAG,valAG);

    return AG;
end

function buildSparseAG_A3DConv(k::Float64,X::Array{Float64,1},Y::Array{Float64,1},
                          Z::Array{Float64,1},
                          fastconv::FastM3D, n::Int64 ,m::Int64,l::Int64; method::String = "normal")
# function that build the sparsigying preconditioner

    # Quick fix: TODO change this!

    Entries = entriesSparseG3D(k,X,Y,Z,fastconv, n ,m, l);
    Ind = reshape(collect(1:n*m*l),n,m,l);

    if method == "normal"
      (Indices, Values) = entriesSparseA3D(k,X,Y,Z,fastconv, n ,m,l);
    elseif method == "randomized"
      println("Not implemented yet")
      (Indices, Values) = entriesSparseARand(k,X,Y,D0, n ,m);
    end

    (rowA, colA, valA) = createIndices(Ind[2:end-1,2:end-1,2:end-1][:],
                                    Indices[1][:], Values[1][:]);

    # building the indices, columns and rows for the interior
    (rowAG, colAG, valAG) = createIndices(Ind[2:end-1,2:end-1,2:end-1][:],
                                    Indices[1][:], (Values[1]*Entries[1])[ :]);

     # for  x = xmin,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,2:end-1][:],
                                    Indices[2][:], Values[2][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmin,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,2:end-1][:],
                                    Indices[2][:], (Values[2]*Entries[2])[ :]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

     # for  x = xmax,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,2:end-1][:],
                                    Indices[3][:], Values[3][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  x = xmax,  y = anything z = anything
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,2:end-1][:],
                                    Indices[3][:], (Values[3]*Entries[3])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    # for  y = ymin,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,2:end-1][:],
                                    Indices[4][:], Values[4][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  y = ymin,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,2:end-1][:],
                                    Indices[4][:], (Values[4]*Entries[4])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    # for  y = ymax,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,2:end-1][:],
                                    Indices[5][:], Values[5][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  y = ymax,  x = anything z = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,2:end-1][:],
                                    Indices[5][:], (Values[5]*Entries[5])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    # for  z = zmin,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,1][:],
                                    Indices[6][:], Values[6][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  z = zmin,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,1][:],
                                    Indices[6][:], (Values[6]*Entries[6])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    # for  z = zmax,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,end][:],
                                    Indices[7][:], Values[7][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    # for  z = zmax,  x = anything y = anything
    (Row, Col, Val) = createIndices(Ind[2:end-1,2:end-1,end][:],
                                    Indices[7][:], (Values[7]*Entries[7])[ :]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

# we need to incorporate the vertices

    # for  x = xmin,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[1,1,2:end-1][:],
                                    Indices[8][:], Values[8][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmin,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[1,1,2:end-1][:],
                                    Indices[8][:], (Values[8]*Entries[8])[ :]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

    # for  x = xmax,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[end,1,2:end-1][:],
                                    Indices[9][:], Values[9][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmax,  y = ymin z = anything
    (Row, Col, Val) = createIndices(Ind[end,1,2:end-1][:],
                                    Indices[9][:], (Values[9]*Entries[9])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

    # for  x = xmin,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,end,2:end-1][:],
                                    Indices[10][:], Values[10][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmin,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,end,2:end-1][:],
                                    Indices[10][:], (Values[10]*Entries[10])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

    # for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[end,end,2:end-1][:],
                                    Indices[11][:], Values[11][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[end,end,2:end-1][:],
                                    Indices[11][:], (Values[11]*Entries[11])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,1][:],
                                    Indices[12][:], Values[12][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmax,  y = ymax z = anything
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,1][:],
                                    Indices[12][:], (Values[12]*Entries[12])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# for  x = xmax,  y = any z = zmin
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,1][:],
                                    Indices[13][:], Values[13][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);
    # for  x = xmax,  y = any z = zmin
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,1][:],
                                    Indices[13][:], (Values[13]*Entries[13])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# for  x = xmin,  y = any z = zmax
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,end][:],
                                    Indices[14][:], Values[14][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);
    # for  x = xmin,  y = any z = zmax
    (Row, Col, Val) = createIndices(Ind[1,2:end-1,end][:],
                                    Indices[14][:], (Values[14]*Entries[14])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# for  x = xmax,  y = any  z = zmax
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,end][:],
                                    Indices[15][:], Values[15][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = xmax,  y = any  z = zmax
    (Row, Col, Val) = createIndices(Ind[end,2:end-1,end][:],
                                    Indices[15][:], (Values[15]*Entries[15])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# for  x = any,  y = tmin  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,1][:],
                                    Indices[16][:], Values[16][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    # for  x = any,  y = tmin  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,1][:],
                                    Indices[16][:], (Values[16]*Entries[16])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex10 = round(Integer, changeInd3D(nHalf,m,1,n,m,l) + subStencil3D(2,1,3,Ind_relative));
# for  x = any,  y = ymax  z = zmin
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,1][:],
                                    Indices[17][:], (Values[17]*Entries[17])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex11 = round(Integer, changeInd3D(nHalf,1,l,n,m,l) + subStencil3D(2,3,1,Ind_relative));
# for  x = any,  y = ymin  z = zmax
    (Row, Col, Val) = createIndices(Ind[2:end-1,1,end][:],
                                    Indices[18][:], (Values[18]*Entries[18])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

# indvertex12 = round(Integer, changeInd3D(nHalf,m,l,n,m,l) + subStencil3D(2,1,1,Ind_relative));
# for  x = any,  y = ymax  z = zmax
    (Row, Col, Val) = createIndices(Ind[2:end-1,end,end][:],
                                    Indices[19][:], (Values[19]*Entries[19])[:]);

    rowAG = vcat(rowAG,Row);
    colAG = vcat(colAG,Col);
    valAG = vcat(valAG,Val);

#################################################################################
############## Now we incorporate the corners  ##################################
#################################################################################

# indcorner1 = round(Integer, changeInd3D(1,1,1,n,m,l) + Ind_relative[2:3,2:3,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[1,1,1],
                                    Indices[20][:], (Values[20]*Entries[20])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner2 = round(Integer, changeInd3D(n,1,1,n,m,l) + Ind_relative[1:2,2:3,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[end,1,1],
                                    Indices[21][:], (Values[21]*Entries[21])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner3 = round(Integer, changeInd3D(1,m,1,n,m,l) + Ind_relative[2:3,1:2,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[1,end,1],
                                    Indices[22][:], (Values[22]*Entries[22])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner4 = round(Integer, changeInd3D(n,m,1,n,m,l) + Ind_relative[1:2,1:2,2:3][:]);
    (Row, Col, Val) = createIndices(Ind[end,end,1],
                                    Indices[23][:], (Values[23]*Entries[23])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner5 = round(Integer, changeInd3D(1,1,l,n,m,l) + Ind_relative[2:3,2:3,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[1,1,end],
                                    Indices[24][:], (Values[24]*Entries[24])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner6 = round(Integer, changeInd3D(n,1,l,n,m,l) + Ind_relative[1:2,2:3,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[end,1,end],
                                    Indices[25][:], (Values[25]*Entries[25])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner7 = round(Integer, changeInd3D(1,m,l,n,m,l) + Ind_relative[2:3,1:2,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[1,end,end],
                                    Indices[26][:], (Values[26]*Entries[26])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)
# indcorner8 = round(Integer, changeInd3D(n,m,l,n,m,l) + Ind_relative[1:2,1:2,1:2][:]);
    (Row, Col, Val) = createIndices(Ind[end,end,end],
                                    Indices[27][:], (Values[27]*Entries[27])[:]);

    rowAG = vcat(rowAG,Row)
    colAG = vcat(colAG,Col)
    valAG = vcat(valAG,Val)

    AG = sparse(rowAG,colAG,valAG);

    return AG;
end

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
module DungAnalyse

using DungBase, ProgressMeter, DelimitedFiles, ProgressMeter
# import Base.Threads: @spawn, @threads

export main, homing, searching, searchcenter, turningpoint

include("load_from_csv.jl")
include("ffmpeg.jl")
include("calibrate.jl")
include("common.jl")

function temp2pixel(coffeesource, temporal2pixel, k, p::POI{C, Temporal{V, I}}) where {C <: Calibration, V <: AbstractTimeLine, I <: Instantaneous}
    xyt = readdlm(joinpath(coffeesource, "pixel", temporal2pixel[Pair(k, p)]))
    Instantaneous(xyt)
end
function temp2pixel(coffeesource, temporal2pixel, k, p::POI{C, Temporal{V, P}}) where {C <: Calibration, V <: AbstractTimeLine, P <: Prolonged}
    xyt = readdlm(joinpath(coffeesource, "pixel", temporal2pixel[Pair(k, p)]))
    Prolonged(xyt)
end

function build_calibrations(coffeesource, temporal2pixel, data)
    cs = unique(p.calib for (_,v) in data for r in v.runs for (_,p) in r.data)
    n = length(cs)
    p = Progress(n, 1, "calibrating")
    Dict((next!(p); c => build_calibration(coffeesource, hash(c), c)) for c in cs)
end
# function build_calibrations(coffeesource, temporal2pixel, data)
#     cs = unique(p.calib for (_,v) in data for r in v.runs for (_,p) in r.data)
#     Dict(c => build_calibration(coffeesource, hash(c), c) for c in unique(p.calib for (_,v) in data for r in v.runs for (_,p) in r.data))
# end

# raw2cm(data, nameerror, coffeesource, temporal2pixel) = Dict(experimentid => Experiment([Run(common(Run(Dict(poitype => calibrate(nameerror[p.calib].filename, temp2pixel(coffeesource, temporal2pixel, poitype, p)) for (poitype, p) in r.data), r.metadata)), r.metadata) for r in v.runs], v.description) for (experimentid, v) in data)
function raw2cm(data, nameerror, coffeesource, temporal2pixel)
    n = 0
    for v in values(data)
        for r in v.runs
            for _ in r.data
                n += 1
            end
        end
    end
    ph = Progress(n)
    trackdata = Dict{String, Experiment}()
    for (k, v) in data
        runs = Vector{Run}(undef, length(v.runs))
        for (i, r) in enumerate(v.runs)
            pois = Dict{Symbol, Any}()
            for (poitype, p) in r.data
                pois[poitype] = calibrate(nameerror[p.calib].filename, temp2pixel(coffeesource, temporal2pixel, poitype, p))
                next!(ph)
            end
            runs[i] = Run(common(Run(pois, r.metadata)), r.metadata)
        end
        trackdata[k] = Experiment(runs, v.description)
    end
    return trackdata
end

function main(coffeesource)
    temporal2pixel, data = loaddeomcsv(coffeesource)
    nameerror = build_calibrations(coffeesource, temporal2pixel, data)
    raw2cm(data, nameerror, coffeesource, temporal2pixel)
end

end # module

# TODO
# change terminology to nest and imaginary nest
# make sure all the terms (like direction etc) are consistent
# define where the binary stops are for the pipeline
# consider how the code works on the experiment[run] versus just run[]
# check all runs
# filter away bad ones
# adjust turning point algorithm
# consider saving the calibration images in a tmp folder instead
# have a thorough walk through the code to fine tune stuff
# consider that calibrations might be better connected to the videos than to the temporals...
# create the regression plot for marie


text:
# fdf(x, α, d) = f(x + α * d), ∇f(x + α * d)
"""
    DescentMethods.strong_backtracking!(fdf, x₀, d, [y₀, g₀]; [α, αmax, β, σ])

Find `α` such that `f(x₀ + αd) ≤ f(x₀) + β α d⋅∇f(x)` and `|d⋅∇f(x₀ + αd)| ≤ σ|d⋅∇f(x₀)|`
    (strong Wolfe conditions) using the backtracking line search with cubic interpolation
    and Hager-Zhang approximation when steps are very small. Optionally, `y₀ = f(x₀)` and
    `g₀ = ∇f(x₀)` may be provided to avoid recalculation.

# Arguments
- `fdf`: a function `fdf(x, α, d)` returning a tuple `f(x + α * d), ∇f(x + α * d)` where `f`
    is the minimized function
- `x₀::AbstractVector`: the initial point
- `d::AbstractVector`: the search direction. Must be a descent direction, i.e. `d⋅∇f(x₀) < 0`
- `y₀`: (optional) the value of `f(x₀)` if it is known beforehand
- `g₀`: (optional) the value of `∇f(x₀)` if it is known beforehand

# Keywords
- `α=1`: the initial value of `α`
- `αmax=Inf`: the maximum allowed value of α
- `β=1e-4`: the coefficient in Wolfe conditions
- `σ=0.5`: the coefficient in Wolfe conditions
"""
function strong_backtracking!(
    fdf::F, x0, d, y0::T, grad0;
    α = one(y0),
    αmax = convert(T, Inf),
    β = convert(T, 1e-4),
    σ = convert(T, 0.5)
) where {F,T}
    α = min(α, αmax / 2)
    α_prev = zero(α)
    y_prev = y0
    ylo = y_prev
    yhi = convert(T, NaN)
    αlo = αhi = convert(T, NaN)
    g0 = dot(grad0, d)
    @assert g0 < 0 "derivative is non-negative: $g0"

    glo = g_prev = g0
    ghi = convert(typeof(g0), NaN)
    # if g0 > 0
    #     g0 = -g0
    #     rmul!(d, -1)
    # end
    mag = max(abs(y0), -α * g0)
    ϵ = sqrt(eps(one(y0))) * mag

    wolfe1 = β * g0
    wolfe2 = -σ * g0

    # bracketing phase
    min_factor = 17/16 # min. factor to expand bracketing interval

    nbracket_max = 200
    for nbracket in 1:nbracket_max
        y, grad = fdf(x0, α, d)
        g = dot(grad, d)
        Δyp = (g + g0) * α / 2 # parabolic approximation
        if abs(Δyp) < ϵ
            Δy = Δyp
        else
            Δy = y - y0
        end
        if Δy > α * wolfe1 || y >= y_prev + ϵ # x >= NaN is always false
            αlo, αhi, ylo, yhi, glo, ghi = α_prev, α, y_prev, y, g_prev, g
            break
        end
        if abs(g) <= wolfe2
            return α
        elseif g >= 0
            αlo, αhi, ylo, yhi, glo, ghi = α_prev, α, y_prev, y, g_prev, g
            break
        end

        nbracket == nbracket_max && error("Failed to find bracketing")
        # cubic interpolation (Nocedal & Wright 2nd ed., p.59)
        Δα = α - α_prev

        # return 0 signaling that bracketing failed
        if iszero(Δα)
            return Δα
        end
        d1 = g_prev + g - 3 * (y - y_prev) / Δα
        det = d1^2 - g * g_prev
        if det < 0
            αnew = min(2, sqrt(αmax / α)) * α
        else
            d2 = sqrt(det)
            αnew = α - Δα * (g + d2 - d1) / (g - g_prev + 2 * d2)
        end
        if α < αnew < αmax && min_factor * α < αnew
            α_prev, α = α, αnew
        else
            α_prev, α = α, α * min(2, sqrt(αmax / α))
        end
        y_prev, g_prev = y, g
    end

    # zoom phase
    # αlo is lower bound, αhi is higher
    small_α = sqrt(eps(one(αlo)))
    for nzoom in 1:200
        # cubic interpolation (Nocedal & Wright 2nd ed., p.59)
        Δα = αhi - αlo
        if Δα < eps(αhi) * 32
            @warn "Step too small; returning"
            return αlo + Δα / 2
        end
        d1 = ghi + glo - 3 * (yhi - ylo) / Δα
        det = d1^2 - ghi * glo
        if det < 0
            α = αlo + Δα / 2
        else
            d2 = sqrt(det)
            α = αlo - Δα * (ghi + d2 - d1) / (ghi - glo + 2 * d2)
            # ensure that new α is in (αlo; αhi) and not too close to bounds
            if !(α - αlo > small_α * αlo && αhi - α > small_α * αhi)
                α = αlo + Δα / 2
            end
        end
        y, grad = fdf(x0, α, d)
        g = grad ⋅ d
        Δyp = (g + g0) * α / 2 # parabolic approximation
        if abs(Δyp) < ϵ
            Δy = Δyp
        else
            Δy = y - y0
        end
        if Δy > α * wolfe1 || y >= ylo + ϵ
            αhi, yhi, ghi = α, y, g
        else
            abs(g) <= wolfe2 && return α
            if g > 0
                αhi, yhi, ghi = α, y, g
            else
                αlo, ylo, glo = α, y, g
            end
        end
    end
end

function strong_backtracking!(
    fdf::F, x0::AbstractVector, d::AbstractVector{T};
    α = one(T),
    αmax = convert(T, Inf),
    β = convert(T, 1e-4),
    σ = convert(T, 0.5)
) where {F,T}
    y0, grad0 = fdf(x0, zero(T), d)
    return strong_backtracking!(fdf, x0, d, y0, grad0, α = α, αmax = αmax, β = β, σ = σ)
end


text:
using GridInterpolations
using LocalFunctionApproximation
using StaticArrays
using JLD2, FileIO
using Random
using ContinuumWorld
using CMSSPs

# Read parameters from file
# params_fn = ARGS[1]
# inhor_fn = ARGS[2]
# outhor_fn = ARGS[3]

rng = MersenneTwister(10)

params_fn = "grids-continuum-params1.toml"
params = continuum_parse_params(params_fn)

# Create CWorld actions
grid_cworld_actions = [Vec2(params.move_amt, 0.0), Vec2(-params.move_amt, 0.0),
                       Vec2(0.0, params.move_amt), Vec2(0.0, -params.move_amt), Vec2(0.0, 0.0)]

# Use one cworld for all of them
cworld = CWorld(xlim = (0.0, 1.0), ylim = (0.0, 1.0),
                reward_regions = [], rewards = [], terminal = [],
                stdev = params.move_std, actions = grid_cworld_actions, discount = 1.0)
cworlds = Dict(1=>cworld, 2=>cworld, 3=>cworld, 4=>cworld)
params.cworlds = cworlds

# Only need to find policy for one mode, since shared
modal_mdp = GridsContinuumMDPType(1, params, grid_cworld_actions, params.beta, params.horizon_limit)

# Create grid interpolator
xy_spacing = polyspace_symmetric(1.0, params.vals_along_axis, 3)
xy_grid = RectangleGrid(xy_spacing, xy_spacing)
lfa = LocalGIFunctionApproximator(xy_grid)

# Compute and update the terminal cost
compute_terminalcost_localapprox!(modal_mdp, lfa)

# Compute policy with localApproxVI
modal_policy = finite_horizon_VI_localapprox!(modal_mdp, lfa, true,
                                              params.num_generative_samples, rng)

# Compute min_cost_per_horizon
compute_min_value_per_horizon_localapprox!(modal_policy)

# Save horizon policy to file
inhor_fn = "grids-continuum-params1-inhor.jld2"
outhor_fn = "grids-continuum-params1-outhor.jld2"
save_modal_horizon_policy_localapprox(modal_policy, inhor_fn,outhor_fn, modal_mdp)


text:
using Test
using SafeTestsets

@safetestset "Geometry" begin include("Geometry/runtests.jl") end

@safetestset "Integration" begin include("Integration/runtests.jl") end

@safetestset "Mesh" begin include("Mesh/runtests.jl") end

# @safetestset "FEM" begin include("FEM/runtests.jl") end

@safetestset "Nystrom" begin include("Nystrom/runtests.jl") end

@safetestset "IO" begin include("IO/runtests.jl") end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open OpToStr.jl")
TabLevel=TabLevel*"    "

###Pretty print for the operatos
function OpPrintln(Op;tex=false,bl=false)
    println(OpToStr(Op;tex=tex,bl=bl))
end
function OpPrint(Op;tex=false,bl=false)
    print(OpToStr(Op;tex=tex,bl=bl))
end

function OpToStr(Num::Number;tex::Bool=false,bl::Bool=false)
    return "$Num"
end
function OpToStr(OP::MajoranaSp;tex::Bool=false,bl::Bool=false)
    n=0
    Str=""
    for B in OP.SpList
        if n==0
            Str *= OpToStr(B,tex=tex)
        else
            if tex
                Str *= " "
            elseif bl
                Str *= "\n"
            else
                Str *= "  "
            end
            Str *= OpToStr(B,tex=tex,FirstTerm=false)
        end
        n+=1
    end
    return Str
end

function OpRealToStr(scale,KetState,NumOps,FirstTerm)
    Rscale=real(scale)
    Str=""
    if Rscale == 1
        ###do nothing except if the list is
        if FirstTerm
            if NumOps==0 && !KetState
                Str*="1"
            end
        else
            Str*="+"
        end
    elseif Rscale ==-1
        if NumOps==0 && !KetState
            Str*="-1"
        else
            Str*="-"
        end
    else
        if FirstTerm ;
            Str*=string(Rscale)*" "
        elseif Rscale >=0
            Str*="+"*string(Rscale)*" "
        else
            Str*=string(Rscale)*" "
        end
    end
    return Str
end

function OpImagToStr(scale,tex,FirstTerm)
    ImScale=imag(scale)
    Str=""
    if ImScale >= 0
        if FirstTerm
            if tex
                Str*=string(ImScale)*"\\i "
            else
                Str*=string(ImScale)*"i "
            end
        else
            if tex
                Str*="+"*string(ImScale)*"\\i "
            else
                    Str*="+"*string(ImScale)*"i "
            end
        end
    else
        if tex
            Str*=string(ImScale)*"\\i "
        else
            Str*=string(ImScale)*"i "
        end
    end
end

function ImStr(FirstTerm,tex)
    if FirstTerm
        if tex
            Str="\\i "
        else
            Str="i "
            end
    else
        if tex
            Str="+\\i "
        else
            Str="+i "
        end
    end
    return Str
end

function ScaleToStr(scale,tex::Bool,KetState::Bool,NumOps::Int,FirstTerm::Bool)
    Str=""
    if isreal(scale)
        Str*=OpRealToStr(scale,KetState,NumOps,FirstTerm)
    elseif scale == im
        Str*=ImStr(FirstTerm,tex)
    elseif scale ==-im
        if tex
            Str*="-\\i "
        else
            Str*="-i "
        end
    elseif isreal(im*scale)
        Str*=OpImagToStr(scale,tex,FirstTerm)
    else
        if !FirstTerm ; Str*="+" end
        Str*="("*string(scale)*") "
    end
    return Str
end

###Prints the Operator list
function OpToStr(OP::MajoranaOp;tex::Bool=false,FirstTerm=true,bl::Bool=true)
    scale=OP.scale
    NumOps=length(OP.OpList)
    Str=ScaleToStr(scale,tex,OP.KetState,NumOps,FirstTerm)
    for j in 1:NumOps
        Str*=MajToStr(OP.OpList[j],tex=tex)
        if j!=NumOps
            Str*=" "
        end
    end
    if OP.KetState
        Str*=KetStr(tex)
    end
    return Str
end

function KetStr(tex::Bool=false)
    if tex
        return "\\left|0\\right\\rangle"
    else
        return "|0>"
    end
end

function MajToStr(Elem::MajElem;tex::Bool=false)
    A = Elem.Type
    if !tex
        Str=string(A)
    else
        if A=="bx"
            Str="b^x"
        elseif A=="by"
            Str="b^y"
        elseif A=="bz"
            Str="b^z"
        elseif A=="c"
            Str="c"
        else
            error("Unknown symbol: $A")
        end
    end
    if Elem.Index == ""
        return Str
    else
        if tex
            return Str*"_{"*string(Elem.Index)*"}"
        else
            return Str*"_"*string(Elem.Index)
        end
    end
end

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close OpToStr.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "LibCURL"
version = v"7.61.0"

# Collection of sources required to build LibCURL
sources = [
    "https://curl.haxx.se/download/curl-7.61.0.tar.gz" =>
    "64141f0db4945268a21b490d58806b97c615d3d0c75bf8c335bbe0efd13b45b5",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/curl-7.61.0
./configure --prefix=$prefix --host=$target --with-mbedtls --without-ssl --disable-manual
if [[ $target == *-w64-mingw32 ]]; then
    LDFLAGS="-L$prefix/bin"
elif [[ $target == x86_64-apple-darwin14 ]]; then
    LDFLAGS="-L$prefix/lib -Wl,-rpath,$prefix/lib"
else
    LDFLAGS="-L$prefix/lib -Wl,-rpath-link,$prefix/lib"
fi
make -j${nproc} LDFLAGS="$LDFLAGS"
make install-exec

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:powerpc64le, :glibc),
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl, :eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libcurl", :libcurl),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.2/build_Zlib.v1.2.11.jl",
    "https://github.com/JuliaWeb/MbedTLSBuilder/releases/download/v0.11/build_MbedTLS.v1.0.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: fchorney/LibCURLBuilder

max_issues_repo_name: fchorney/LibCURLBuilder

max_forks_repo_name: fchorney/LibCURLBuilder


text:
# ------------------------------------------------------------------

mutable struct node
  level::Int
  model::JuMP.Model
end

## Checks if solution is binary
function isBinary(model::JuMP.Model, binaryIndices::Vector{Int64})
  isZero = model.colVal[binaryIndices] .== 0
  isOne  = model.colVal[binaryIndices] .== 1
  if sum(isZero + isOne) == length(binaryIndices)
    return true
  end
  return false
end

## Checks if model is max: if it is, converts to min
function convertSense!(m::JuMP.Model)
  if m.objSense == :Max
    m.objSense = :Min
    m.obj = -m.obj
  else
    m.objSense = :Max
    m.obj = -m.obj
  end
end

## Receives node and creates two children by setting a variable to 0 and 1 respectively
function branch(currentNode::node, binaryIndices::Vector{Int64})

  distance = abs(currentNode.model.colVal[binaryIndices] - 0.5)
  indFrac = indmin(distance)
  indToSet = binaryIndices[indFrac]

  leftModel = deepcopy(currentNode.model)
  leftModel.colUpper[indToSet] = 0
  leftModel.colLower[indToSet] = 0

  rightModel = deepcopy(currentNode.model)
  rightModel.colUpper[indToSet] = 1
  rightModel.colLower[indToSet] = 1

  leftChild = node(currentNode.level+1, leftModel)
  rightChild = node(currentNode.level+1, rightModel)

  return leftChild, rightChild
end

function obtainBoundList(nodeList::Vector{node})
  boundList = Array{Float64}(length(nodes))
  for i = 1 : length(nodes)
    boundList[i] = nodes[i].model.objVal
  end

  return boundList
end

## Receives a mixed binary linear JuMP model
function solveMIP(m::JuMP.Model)

  tic()

  # Check if model is max; if it is, converts to min
  flagConverted = 0
  if m.objSense == :Max
    convertSense!(m)
    flagConverted = 1
  end

  # Best bounds: start as Inf
  bestBound = -1e200
  bestVal = 1e200

  # Create vector of indices of the binary variables
  binaryIndices = find(m.colCat .== :Bin)
  binarySolutions = 0

  # Solve linear relaxation
  status = solve(m, relaxation=true)
  if status == :Optimal && isBinary(m, binaryIndices)
    # Solution of the relaxed problem is binary: optimal solution
    nodes = Vector{node}(0)
    bestBound = m.objVal
    bestVal = m.objVal
    binarySolutions = 1
  else
    # Create branch and bound tree
    nodes = Vector{node}(1)
    nodes[1] = node(0, m) # root
    lastNodeLevel = 0
  end

  iter = 0
  flagOpt = 0
  tol = 0.01 # Tolerance (%)
  time0 = time_ns()
  while !isempty(nodes) && abs((bestVal - bestBound)/bestVal) > tol && (time_ns()-time0)/1e9 < 600
    if iter == 0
      iter = 1
      continue
    end
    # Check node lower bound. If greater than current best UB, prune by limit
    if nodes[1].model.objVal <= bestVal
      status = solve(nodes[1].model)
      if status == :Optimal
        bestBound = minimum(obtainBoundList(nodes))
        if isBinary(nodes[1].model, binaryIndices)
          # Relaxed solution is binary: optimal solution -- don't branch
          if nodes[1].model.objVal < bestVal
            bestVal = nodes[1].model.objVal
            m.colVal = copy(nodes[1].model.colVal)
            flagOpt = 1
            binarySolutions+=1
          end
        elseif nodes[1].model.objVal <= bestVal
          # Relaxed solution is not binary and should not be pruned by limit -- branch
          (leftChild, rightChild) = branch(nodes[1], binaryIndices)
          push!(nodes, leftChild)
          push!(nodes, rightChild)
        end
      end
    end
    lastNodeLevel = nodes[1].level
    deleteat!(nodes, 1)
    iter+=1

    if iter == 1 || iter%10 == 0
      println("UB: $bestVal")
      println("LB: $bestBound")
    end
  end

  if flagConverted == 1
    convertSense!(m)
    m.objVal = -bestVal
    m.objBound = -bestBound
  else
    m.objVal = bestVal
    m.objBound = bestBound
  end

  # Outputs
  m.ext[:status] = status
  if flagOpt == 1
    m.ext[:status] = :Optimal
  end
  m.ext[:nodes] = iter
  m.ext[:solutions] = binarySolutions
  t = toc()
  m.ext[:time] = t

  return m.ext[:status]
end


text:
Base.@kwdef struct Config
    pkg_server::String
    before_script::String = _DEFAULT_BEFORE_SCRIPT
    script::String        = _DEFAULT_SCRIPT
    after_script::String  = _DEFAULT_AFTER_SCRIPT
end

Base.@kwdef struct TemporaryEnvironment
    env_with_project::Dict{String, String}
    env_without_project::Dict{String, String}
end

Base.@kwdef struct Commands
    cmd_1_before_script::Cmd
    cmd_2_script::Cmd
    cmd_3_after_script::Cmd
end

Base.@kwdef struct ProjectFileOrManifestFile
    contents::String
    name::String
    before_script::Union{String, Nothing} = nothing
    script::Union{String, Nothing} = nothing
    after_script::Union{String, Nothing} = nothing
end

function Project(contents::AbstractString;
                 before_script = nothing,
                 script = nothing,
                 after_script = nothing)
    name = "Project.toml"
    result = ProjectFileOrManifestFile(;
        contents,
        name,
        before_script,
        script,
        after_script,
    )
    return result
end

function Manifest(contents::AbstractString;
                 before_script = nothing,
                 script = nothing,
                 after_script = nothing)
    name = "Manifest.toml"
    result = ProjectFileOrManifestFile(;
        contents,
        name,
        before_script,
        script,
        after_script,
    )
    return result
end

FileName(filename::AbstractString)::String = read(filename, String)::String
FileContents(contents::AbstractString)::String = convert(String, contents)::String


text:
module HADSGE
using Calculus,Base.Threads
import Base:names,length,getindex,display,find,parse,spzeros,kron,*,promote_type,setindex!,Future
import SparseGrids:NGrid,ndgrid,@threadsfixed
using SparseGrids
include("utils.jl")
include("parse.jl")
include("calculus.jl")
include("parsemodel.jl")
include("model.jl")
include("initmodel.jl")
include("solve.jl")
include("aggregate.jl")
include("simulation.jl")
include("adapt.jl")
# include("display.jl")

export Model, solve, ∫, updateA

end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE