### Init some initial parameters necessary to create the plan ###

function initParams(x::Matrix{T}, N::NTuple{D,Int}, dims::Union{Integer,UnitRange{Int64}}=1:D; 
                    kargs...) where {D,T}
  # convert dims to a unit range
  dims_ = (typeof(dims) <: Integer) ? (dims:dims) : dims

  params = NFFTParams{T,D}(; kargs...)
  m, σ, reltol = accuracyParams(; kargs...)
  params.m = m
  params.σ = σ
  params.reltol = reltol

  # Taken from NFFT3
  m2K = [1, 3, 7, 9, 14, 17, 20, 23, 24]
  K = m2K[min(m+1,length(m2K))]
  params.LUTSize = 2^(K) * (m) # ensure that LUTSize is dividable by (m)

  if length(dims_) != size(x,1)
      throw(ArgumentError("Nodes x have dimension $(size(x,1)) != $(length(dims_))"))
  end

  if any(isodd.(N[dims]))
    throw(ArgumentError("N = $N needs to consist of even integers along dims = $(dims)!"))
  end

  doTrafo = ntuple(d->d ∈ dims_, D)

  n = ntuple(d -> doTrafo[d] ?
                      (ceil(Int,params.σ*N[d])÷2)*2 : # ensure that n is an even integer
                        N[d], D)

  params.σ = n[dims_[1]] / N[dims_[1]]

  #params.blockSize = ntuple(d-> n[d] , D) # just one block
  if haskey(kargs, :blockSize)
    params.blockSize = kargs[:blockSize]
  else
    params.blockSize = ntuple(d-> _blockSize(n,d) , D)
  end

  M = size(x, 2)

  # calculate output size
  NOut = Int[]
  Mtaken = false
  for d=1:D
    if !doTrafo[d]
      push!(NOut, N[d])
    elseif !Mtaken
      push!(NOut, M)
      Mtaken = true
    end
  end
  # Sort nodes in lexicographic way
  if params.sortNodes
      x .= sortslices(x, dims=2)
  end
  return params, N, Tuple(NOut), M, n, dims_
end


function _blockSize(n::NTuple{1,Int}, d)
  return min(1024, n[d])
end

function _blockSize(n::NTuple{2,Int}, d)
  return min(64, n[d])
end

function _blockSize(n::NTuple{D,Int}, d) where {D}
  if d == 1
    return min(16, n[d])
  elseif d == 2
    return min(16, n[d])
  elseif d == 3
    return min(16, n[d])
  else
    return 1
  end
end

### Precomputation of the B matrix ###

function precomputeB(win, x, N::NTuple{D,Int}, n::NTuple{D,Int}, m, M, σ, K, T) where D
  I = Array{Int64,2}(undef, (2*m)^D, M)
  β = (2*m)^D
  J = [β*k+1 for k=0:M]
  V = Array{T,2}(undef,(2*m)^D, M)
  mProd = ntuple(d-> (d==1) ? 1 : (2*m)^(d-1), D)
  nProd = ntuple(d-> (d==1) ? 1 : prod(n[1:(d-1)]), D)
  L = Val(2*m)
  scale = Int(K/(m))

  @cthreads for k in 1:M
    _precomputeB(win, x, N, n, m, M, σ, scale, I, J, V, mProd, nProd, L, k, K)
  end

  S = SparseMatrixCSC(prod(n), M, J, vec(I), vec(V))
  return S
end

@inline @generated function _precomputeB(win, x::AbstractMatrix{T}, N::NTuple{D,Int}, n::NTuple{D,Int}, m, M,
                     σ, scale, I, J, V, mProd, nProd, L::Val{Z}, k, LUTSize) where {T, D, Z}
  quote

    @nexprs $(D) d -> ((tmpIdx_d, tmpWin_d) = precomputeOneNode(win, x, n, m, σ, scale, k, d, L, LUTSize) )

    @nexprs 1 d -> κ_{$D} = 1 # This is a hack, I actually want to write κ_$D = 1
    @nexprs 1 d -> ζ_{$D} = 1
    @nexprs 1 d -> prodWin_{$D} = one(T)
    @nloops $D l d -> 1:$Z d->begin
        # preexpr
        prodWin_{d-1} = prodWin_d * tmpWin_d[l_d]
        κ_{d-1} = κ_d + (l_d-1) * mProd[d]
        ζ_{d-1} = ζ_d + (tmpIdx_d[l_d]-1) * nProd[d]
    end begin
        # bodyexpr
        I[κ_0,k] = ζ_0
        V[κ_0,k] = prodWin_0
    end
    return
  end
end

### precomputation of the window and the indices required during convolution ###

@generated function precomputeOneNode(win::Function, x::AbstractMatrix{T}, n::NTuple{D,Int}, m,
  σ, scale, k, d, L::Val{Z}, LUTSize) where {T,D,Z}
  quote
    xscale = x[d,k] * n[d]
    off = floor(Int, xscale) - m + 1
    tmpIdx = @ntuple $(Z) l -> ( rem(l + off + n[d] - 1, n[d]) + 1)

    tmpWin = @ntuple $(Z) l -> (win( (xscale - (l-1) - off)  / n[d], n[d], m, σ) )
    return (tmpIdx, tmpWin)
  end
end

# precompute = LINEAR

@generated function precomputeOneNode(winLin::Array, winPoly::Nothing, x::AbstractMatrix{T}, n::NTuple{D,Int}, m,
  σ, scale, k, d, L::Val{Z}, LUTSize) where {T,D,Z}
  quote
    xscale = x[d,k] * n[d]
    off = floor(Int, xscale) - m + 1
    tmpIdx = @ntuple $(Z) l -> ( rem(l + off + n[d] - 1, n[d]) + 1)

    idx = ((xscale - off)*LUTSize)/(m)
    tmpWin =  shiftedWindowEntries(winLin, idx, scale, d, L)

    return (tmpIdx, tmpWin)
  end
end

# precompute = POLYNOMIAL

@generated function precomputeOneNode(winLin::Array, winPoly::NTuple{Y, NTuple{X,T}}, x::AbstractMatrix{T}, n::NTuple{D,Int}, m,
  σ, scale, k, d, L::Val{Z}, LUTSize) where {Y,X,T,D,Z}
  quote
    xscale = x[d,k] * n[d]
    off = floor(Int, xscale) - m + 1
    tmpIdx = @ntuple $(Z) l -> ( rem(l + off + n[d] - 1, n[d]) + 1)

    idx = (xscale - off - m + T(0.5))
    tmpWin =  shiftedWindowEntries(winPoly, idx, scale, d, L)

    return (tmpIdx, tmpWin)
  end
end

# precompute = LINEAR

function precomputeOneNodeBlocking(winLin, winTensor::Nothing, winPoly::Nothing,
                                    scale, k, d, L, idxInBlock::Matrix)

  y, idx = idxInBlock[d,k]
  tmpWin =  shiftedWindowEntries(winLin, idx, scale, d, L)

  return (y, tmpWin)
end

@generated function shiftedWindowEntries(winLin::Vector, idx, scale, d, L::Val{Z}) where {Z}
  quote
    idxL = floor(Int,idx)
    idxInt = Int(idxL)
    α = ( idx-idxL )

    tmpWin = @ntuple $(Z) l -> begin
      # Uncommented code: This is the version where we pull in l into the abs.
      # We pulled this out of the iteration.
      # idx = abs((xscale - (l-1)  - off)*LUTSize)/(m)

      # The second +1 is because Julia has 1-based indexing
      # The first +1 is part of the index calculation and needs(!)
      # to be part of the abs. The abs is shifting to the positive interval
      # and this +1 matches the `floor` above, which rounds down. In turn
      # for positive and negative indices a different neighbour is calculated
      idxInt1 = abs( idxInt - (l-1)*scale ) +1
      idxInt2 = abs( idxInt - (l-1)*scale +1) +1

      (winLin[idxInt1] + α * (winLin[idxInt2] - winLin[idxInt1]))
    end
    return tmpWin
  end
end

# precompute = POLYNOMIAL

function precomputeOneNodeBlocking(winLin, winTensor::Nothing, winPoly::NTuple{Z, NTuple{X,T}}, scale,
                                   k, d, L, idxInBlock::Matrix) where {T,Z,X}

  y, x = idxInBlock[d,k]
  tmpWin =  shiftedWindowEntries(winPoly, x, scale, d, L)

  return (y, tmpWin)
end


@generated function shiftedWindowEntries(winPoly::NTuple{Z, NTuple{X,T}}, x::T, scale, d, L::Val{Z}) where {T,Z,X}
  quote
    tmpWin = @ntuple $(Z) l -> begin
      evalpoly(x, winPoly[l])
    end
    return tmpWin
  end
end

#= Static Array version. Not faster
function precomputeOneNodeBlocking(winLin, winTensor::Nothing, winPoly::SMatrix{X,Z,T}, scale,
  k, d, L, idxInBlock::Matrix) where {T,Z,X}

  y, x = idxInBlock[d,k]
  tmpWin =  shiftedWindowEntries(winPoly, x, scale, d, L)

  return (y, tmpWin)
end

@generated function shiftedWindowEntries(winPoly::SMatrix{X,Z,T}, x::T, scale, d, L::Val{Z}) where {T,Z,X}
  quote

    x_1 = one(T)
    @nexprs $(X-1) h->(x_{h+1} = x_h * x)

    xx = @ntuple $(X) h -> begin
      x_{h}
    end

    xs = SVector(xx...)

    tmpWin = transpose(winPoly) * xs

    return tmpWin
  end
end
=#

# precompute = TENSOR

function precomputeOneNodeBlocking(winLin, winTensor::Array, winPoly::Nothing, scale,
                                   k, d, L, idxInBlock::Matrix)
  y, idx = idxInBlock[d,k]
  tmpWin =  shiftedWindowEntriesTensor(winTensor, k, d, L)

  return (y, tmpWin)
end


@generated function shiftedWindowEntriesTensor(winTensor, k, d, L::Val{Z}) where {Z}
  quote
    tmpWin = @ntuple $(Z) l -> begin
      winTensor[l, d, k]
    end
    return tmpWin
  end
end


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

## function nfft_precompute_lin_psi in NFFT3
"""
Precompute the look up table for the window function φ.

Remarks:
* Only the positive half is computed
* The window is computed for the interval [0, (m)/n]. The reason for the +2 is
  that we do evaluate the window function outside its interval, since x does not
  necessary match the sampling points
* The window has K+1 entries and during the index calculation we multiply with the
  factor K/(m).
* It is very important that K/(m) is an integer since our index calculation exploits
  this fact. We therefore always use `Int(K/(m))`instead of `K÷(m)` since this gives
  an error while the later variant would silently error.
"""
function precomputeLinInterp(win, m, σ, K, T)
  windowLinInterp = Vector{T}(undef, K+1)

  step = (m) / (K)
  @cthreads for l = 1:(K+1)
      y = ( (l-1) * step )
      windowLinInterp[l] = win(y, 1, m, σ)
  end
  return windowLinInterp
end

function precomputePolyInterp(win, m, σ, T)
  deg = 2*m # 2*m #+ 2   # Certainly depends on Window
  K = 2*m
  NSamples = 2*deg # Sample more densely!!!
  windowPolyInterp = Matrix{T}(undef, deg, K)

  x = range(-0.5, 0.5,length=NSamples)
  V = ones(deg, NSamples)
  for r=2:deg
    V[r,:] .= V[r-1,:] .* x
  end

  for l = 1:K
      y = (-(l-0.5) + m) .+ x
      samples = win.(y, 1, m, σ)
      windowPolyInterp[:,l] .= V' \ samples
  end
  return windowPolyInterp
end

function testPrecomputePoly(win, m, σ, T)
  deg = 2*m + 3
  K = 2*m
  step = 1

  windowPolyInterp = precomputePolyInterp(win, m, σ, T)

  x = -0.394234
  xx = [x^l for l=0:(deg-1)]

  winTrue = zeros(T,K)
  winApprox = zeros(T,K)

  for l = 1:K
    y = (-(l-0.5) * step + m) + x
    @info y
    winTrue[l] = win(y, 1, m, σ)
    winApprox[l] = windowPolyInterp[:,l]' * xx
  end

  return winTrue, winApprox
end


function precomputeWindowHatInvLUT(windowHatInvLUT, win_hat, N, n, m, σ, T)

  for d=1:length(windowHatInvLUT)
      κ = k -> win_hat(k-1-N[d]÷2, n[d], m, σ)
      cheb = ChebyshevInterpolator(κ, 1, N[d], 30)

      windowHatInvLUT[d] = zeros(T, N[d])
      @cthreads for k=1:N[d]
          windowHatInvLUT[d][k] = 1. / cheb(k) #win_hat(k-1-N[d]÷2, n[d], m, σ)
      end
  end
end

function precomputation(x::Union{Matrix{T},Vector{T}}, N::NTuple{D,Int}, n, params) where {T,D}

  m = params.m; σ = params.σ; window=params.window
  LUTSize = params.LUTSize; precompute = params.precompute

  win, win_hat = getWindow(window) # highly type instable. But what should be do
  M = size(x, 2)

  windowHatInvLUT_ = Vector{Vector{T}}(undef, D)
  precomputeWindowHatInvLUT(windowHatInvLUT_, win_hat, N, n, m, σ, T)

  if params.storeDeconvolutionIdx
    windowHatInvLUT = Vector{Vector{T}}(undef, 1)
    windowHatInvLUT[1], deconvolveIdx = precompWindowHatInvLUT(params, N, n, windowHatInvLUT_)
  else
    windowHatInvLUT = windowHatInvLUT_
    deconvolveIdx = Array{Int64,1}(undef, 0)
  end

  if precompute == LINEAR
    windowLinInterp = precomputeLinInterp(win, m, σ, LUTSize, T)
    windowPolyInterp = Matrix{T}(undef, 0, 0)
    B = sparse([],[],T[])
  elseif precompute == POLYNOMIAL
    windowLinInterp = Vector{T}(undef, 0)
    windowPolyInterp = precomputePolyInterp(win, m, σ, T)
    B = sparse([],[],T[])
  elseif precompute == FULL
    windowLinInterp = Vector{T}(undef, 0)
    windowPolyInterp = Matrix{T}(undef, 0, 0)
    B = precomputeB(win, x, N, n, m, M, σ, LUTSize, T)
    #windowLinInterp = precomputeLinInterp(win, windowLinInterp, n, m, σ, LUTSize, T) # These versions are for debugging
    #B = precomputeB(windowLinInterp, x, N, n, m, M, σ, LUTSize, T)
  elseif precompute == TENSOR
    windowLinInterp = Vector{T}(undef, 0)
    windowPolyInterp = Matrix{T}(undef, 0, 0)
    B = sparse([],[],T[])
  else
    windowLinInterp = Vector{T}(undef, 0)
    windowPolyInterp = Matrix{T}(undef, 0, 0)
    B = sparse([],[],T[])
    error("precompute = $precompute not supported by NFFT.jl!")
  end

  return (windowLinInterp, windowPolyInterp, windowHatInvLUT, deconvolveIdx, B)
end

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


# This function is type unstable. why???
function precompWindowHatInvLUT(p::NFFTParams{T}, N, n, windowHatInvLUT_) where {T}

  windowHatInvLUT = zeros(Complex{T}, N)
  deconvIdx = zeros(Int64, N)

  if length(N) == 1
    precompWindowHatInvLUT(p, windowHatInvLUT, deconvIdx, N, n, windowHatInvLUT_, 1)
  else
    @cthreads for o = 1:N[end]
      precompWindowHatInvLUT(p, windowHatInvLUT, deconvIdx, N, n, windowHatInvLUT_, o)
    end
  end
  return vec(windowHatInvLUT), vec(deconvIdx)
end

@generated function precompWindowHatInvLUT(p::NFFTParams{T}, windowHatInvLUT::AbstractArray{Complex{T},D},
           deconvIdx::AbstractArray{Int,D}, N, n, windowHatInvLUT_, o)::Nothing where {D,T}
  quote
    linIdx = LinearIndices(n)

    @nexprs 1 d -> gidx_{$D-1} = rem(o+n[$D] - N[$D]÷2 - 1, n[$D]) + 1
    @nexprs 1 d -> l_{$D-1} = o
    @nloops $(D-2) l d->(1:N[d+1]) d-> begin
        gidx_d = rem(l_d+n[d+1] - N[d+1]÷2 - 1, n[d+1]) + 1
      end begin
      N2 = N[1]÷2
      @inbounds @simd for i = 1:N2
        deconvIdx[i, CartesianIndex(@ntuple $(D-1) l)] =
           linIdx[i-N2+n[1], CartesianIndex(@ntuple $(D-1) gidx)]
        v = windowHatInvLUT_[1][i]
        @nexprs $(D-1) d -> v *= windowHatInvLUT_[d+1][l_d]
        windowHatInvLUT[i, CartesianIndex(@ntuple $(D-1) l)] = v

        deconvIdx[i+N2, CartesianIndex(@ntuple $(D-1) l)] =
           linIdx[i, CartesianIndex(@ntuple $(D-1) gidx)]
        v = windowHatInvLUT_[1][i+N2]
        @nexprs $(D-1) d -> v *= windowHatInvLUT_[d+1][l_d]
        windowHatInvLUT[i+N2, CartesianIndex(@ntuple $(D-1) l)] = v
      end
    end
    return
  end
end

####### block precomputation #########

function precomputeBlocks(x::Matrix{T}, n::NTuple{D,Int}, params, calcBlocks::Bool) where {T,D}

  if calcBlocks
    xShift = copy(x)
    shiftNodes!(xShift)
    blocks, nodesInBlocks, blockOffsets =
        _precomputeBlocks(xShift, n, params.m, params.LUTSize, params.blockSize)

    idxInBlock =  _precomputeIdxInBlock(xShift, n, params.m, params.precompute, params.LUTSize, blockOffsets, nodesInBlocks)
    if params.precompute != TENSOR
      windowTensor = Array{Array{T,3},D}(undef, ntuple(d->0,D))
    else
      #idxInBlock = Array{Matrix{Tuple{Int,Float64}},D}(undef, ntuple(d->0,D))
      windowTensor = _precomputeWindowTensor(xShift, n, params.m, params.σ, nodesInBlocks, params.window)
    end

  else
    blocks = Array{Array{Complex{T},D},D}(undef,ntuple(d->0,D))
    nodesInBlocks = Array{Vector{Int64},D}(undef,ntuple(d->0,D))
    blockOffsets = Array{NTuple{D,Int64},D}(undef,ntuple(d->0,D))
    idxInBlock = Array{Matrix{Tuple{Int,T}},D}(undef, ntuple(d->0,D))
    windowTensor = Array{Array{T,3},D}(undef, ntuple(d->0,D))
  end

  return (blocks, nodesInBlocks, blockOffsets, idxInBlock, windowTensor)
end


function _precomputeBlocks(x::Matrix{T}, n::NTuple{D,Int}, m, LUTSize, blockSize) where {T,D}

  padding = ntuple(d->m, D)
  numBlocks =  ntuple(d-> ceil(Int, n[d]/blockSize[d]), D)
  blockSizePadded = ntuple(d-> blockSize[d] + 2*padding[d], D)
  nodesInBlock = [ Int[] for l in CartesianIndices(numBlocks) ]
  numNodesInBlock = zeros(Int, numBlocks)
  for k=1:size(x,2) # @cthreads
    idx = ntuple(d->unsafe_trunc(Int, x[d,k]*n[d])÷blockSize[d]+1, D)
    numNodesInBlock[idx...] += 1
  end
  @cthreads  for l in CartesianIndices(numBlocks)
    sizehint!(nodesInBlock[l], numNodesInBlock[l])
  end
  for k=1:size(x,2) # @cthreads
    idx = ntuple(d->unsafe_trunc(Int, x[d,k]*n[d])÷blockSize[d]+1, D)
    push!(nodesInBlock[idx...], k)
  end

  blocks = Array{Array{Complex{T},D},D}(undef, numBlocks)
  blockOffsets = Array{NTuple{D,Int64},D}(undef, numBlocks)

  @cthreads for l in CartesianIndices(numBlocks)
    if !isempty(nodesInBlock[l])
      # precompute blocks
      blocks[l] = Array{Complex{T},D}(undef, blockSizePadded)

      # precompute blockOffsets
      blockOffsets[l] = ntuple(d-> (l[d]-1)*blockSize[d]-padding[d]-1, D)
    end
  end

  return blocks, nodesInBlock, blockOffsets
end



function _precomputeIdxInBlock(x::Matrix{T}, n::NTuple{D,Int}, m, precompute, LUTSize, blockOffsets, nodesInBlock) where {T,D}

  numBlocks = size(nodesInBlock)

  idxInBlock = Array{Matrix{Tuple{Int,T}},D}(undef, numBlocks)

  @cthreads for l in CartesianIndices(numBlocks)
    if !isempty(nodesInBlock[l])

      # precompute idxInBlock
      idxInBlock[l] = Matrix{Tuple{Int,T}}(undef, D, length(nodesInBlock[l]))
      @inbounds for (i,k) in enumerate(nodesInBlock[l])
        @inbounds for d=1:D
          xtmp = x[d,k] # this is expensive because of cache misses
          xscale = xtmp * n[d]
          off = unsafe_trunc(Int, xscale) - m + 1
          y = off - blockOffsets[l][d] - 1

          if precompute == LINEAR
            idx = (xscale - off)*(LUTSize÷(m))
          else
            idx = (xscale - off - m + 1 -0.5 )
          end
          idxInBlock[l][d,i] = (y,idx)
        end
      end
    end
  end

  return idxInBlock
end

function _precomputeWindowTensor(x::Matrix{T}, n::NTuple{D,Int}, m, σ, nodesInBlock, window::Symbol) where {T,D}
  win, win_hat = getWindow(window) # highly type instable. But what should be do
  P = precomputePolyInterp(win, m, σ, T)
  winPoly = ntuple(d-> ntuple(g-> P[g,d], size(P,1)), size(P,2))

  return _precomputeWindowTensor(x, n, m, σ, nodesInBlock, winPoly)
end


function _precomputeWindowTensor(x::Matrix{T}, n::NTuple{D,Int}, m, σ, nodesInBlock, winPoly::NTuple{Z, NTuple{X,T}}) where {T,D,Z,X}

  numBlocks = size(nodesInBlock)
  windowTensor = Array{Array{T,3},D}(undef, numBlocks)

  @cthreads for l in CartesianIndices(numBlocks)
    if !isempty(nodesInBlock[l])

      # precompute idxInBlock
      windowTensor[l] = Array{T,3}(undef, 2*m, D, length(nodesInBlock[l]))
      @inbounds for (i,k) in enumerate(nodesInBlock[l])
        @inbounds for d=1:D
          xtmp = x[d,k]  # this is expensive because of cache misses
          xscale = xtmp * n[d]
          off = unsafe_trunc(Int, xscale) - m + 1
          x_ = (xscale - off - m + 1 -0.5 )

          @inbounds @simd for k=1:Z
            windowTensor[l][k,d,i] = evalpoly(x_, winPoly[k])
          end
        end
      end
    end
  end

  return windowTensor
end
