module Numeric.FFTW.Extra.Rank1.Convolution (
   convolveCyclic,
   convolveCyclicOneMany,
   convolveCyclicManyMany,
   convolve,
   convolveShortLong,
   convolveShortLongCut,
   convolveLongToShort,
   convolveLongToShortInterval,
   convolveOneMany,
   convolveManyMany,
   convolveOneManyPad,
   convolveManyManyPad,

   -- export for tests
   arrange,
   arrangeMutable,
   arrangeViaMatrixAdd,
   ) where

import qualified Numeric.FFTW.Extra.Shape as CShape
import Numeric.FFTW.Extra.Utility

import qualified Numeric.FFTW.Rank1 as Rank1
import qualified Numeric.FFTW.Batch as Batch
import qualified Numeric.FFTW.Shape as Spectrum
import qualified Numeric.BLAS.Matrix.RowMajor.Mutable as MutMatrix
import qualified Numeric.BLAS.Matrix.RowMajor.Block as BlockMatrix
import qualified Numeric.BLAS.Matrix.RowMajor as Matrix
import qualified Numeric.BLAS.Vector.Mutable as MutVector
import qualified Numeric.BLAS.Vector.Chunk as Chunk
import qualified Numeric.BLAS.Vector.Slice as VectorSlice
import qualified Numeric.BLAS.Vector as Vector
import qualified Numeric.BLAS.Subobject.Shape as Subshape
import qualified Numeric.BLAS.Subobject.View.Matrix as ViewMatrix
-- import qualified Numeric.BLAS.Subobject.View as View
import qualified Numeric.BLAS.Subobject.Layout as Layout
import qualified Numeric.BLAS.Slice as Slice
import qualified Numeric.BLAS.Scalar as Scalar
import qualified Numeric.Netlib.Class as Class
import Numeric.FFTW.Extra.NumberTheory (ceilingEfficient, ceilingPowerOfTwo)
import Numeric.BLAS.Matrix.RowMajor.Block ((&|||))
import Numeric.BLAS.Vector.Slice ((|+|), (|-|))

import qualified Data.Array.Comfort.Storable.Mutable as MutArray
import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Boxed as BoxedArray
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable (Array)
import Data.Array.Comfort.Shape ((::+)((::+)))

import Data.Bool.HT (if')
import Data.Complex (Complex((:+)))
import Data.Bits (Bits)

import Control.Monad (zipWithM_)

import GHC.Stack (HasCallStack)


{- $setup
>>> :set -XTypeFamilies
>>> import Test.Numeric.FFTW.Extra.Utility
>>> import Numeric.FFTW.Extra.Rank1.Convolution
>>> import Numeric.FFTW.Extra.Utility (padBlock)
>>> import qualified Numeric.FFTW.Extra.Shape as CShape
>>> import qualified Numeric.BLAS.Rank1.Convolution as Naive
>>> import qualified Numeric.BLAS.Matrix.RowMajor as Matrix
>>> import qualified Numeric.BLAS.Vector as Vector
>>> import Numeric.BLAS.Scalar (RealOf)
>>>
>>> import qualified Data.Array.Comfort.Storable.Dim2 as Array2
>>> import qualified Data.Array.Comfort.Storable as Array
>>> import qualified Data.Array.Comfort.Boxed as BoxedArray
>>> import qualified Data.Array.Comfort.Shape as Shape
>>> import Data.Array.Comfort.Shape ((::+)((::+)))
>>> import Data.Complex (Complex)
>>>
>>> import qualified Test.QuickCheck as QC
>>>
>>> type Number_ = NumberType_
>>> type Real_ = RealOf Number_
>>> type Complex_ = Complex Real_
>>>
>>> real_ :: QC.Gen Real_
>>> real_ = genReal
>>>
>>> complex_ :: QC.Gen Complex_
>>> complex_ = genNumber
>>>
>>> number_ :: QC.Gen Number_
>>> number_ = genNumber
-}


newtype Convolve f g h a = Convolve {getConvolve :: f a -> g a -> h a}

{- |
prop> :{
   QC.forAll (genCyclic number_) $ \x ->
   QC.forAll (genVectorForShape number_ $ Array.shape x) $ \y ->
   convolveCyclic x y =~= Naive.convolveCyclic y x
:}
-}
convolveCyclic ::
   (Integral n, Class.Floating a) =>
   Array (Shape.Cyclic n) a ->
   Array (Shape.Cyclic n) a ->
   Array (Shape.Cyclic n) a
convolveCyclic =
   getConvolve $
   Class.switchFloating
      (Convolve convolveCyclicReal)
      (Convolve convolveCyclicReal)
      (Convolve convolveCyclicComplex)
      (Convolve convolveCyclicComplex)

{- |
prop> :{
   QC.forAll (genCyclic number_) $ \x ->
   let width = Array.shape x in
   QC.forAll genVectorShape $ \height ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \ys ->

   convolveCyclicOneMany x ys
   =~=
   Array2.fromRowArray width (fmap (convolveCyclic x) (Array2.toRowArray ys))
:}

prop> :{
   QC.forAll (genCyclic number_) $ \x ->
   let width = Array.shape x in
   QC.forAll genVectorShape $ \height ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \ys ->

   convolveCyclicOneMany x ys
   =~=
   Naive.convolveCyclicOneMany x ys
:}
-}
convolveCyclicOneMany ::
   (Shape.C sh, Integral n, Class.Floating a) =>
   Array     (Shape.Cyclic n) a ->
   Array (sh, Shape.Cyclic n) a ->
   Array (sh, Shape.Cyclic n) a
convolveCyclicOneMany =
   getConvolve $
   Class.switchFloating
      (Convolve convolveCyclicOneManyReal)
      (Convolve convolveCyclicOneManyReal)
      (Convolve convolveCyclicOneManyComplex)
      (Convolve convolveCyclicOneManyComplex)

{- |
prop> :{
   QC.forAll genVectorShape $ \height ->
   QC.forAll genCyclicShape $ \width ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \xs ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \ys ->

   convolveCyclicManyMany xs ys
   =~=
   Array2.fromRowArray width
      (BoxedArray.zipWith convolveCyclic
         (Array2.toRowArray xs) (Array2.toRowArray ys))
:}

prop> :{
   QC.forAll genVectorShape $ \height ->
   QC.forAll genCyclicShape $ \width ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \xs ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \ys ->

   convolveCyclicManyMany xs ys
   =~=
   Naive.convolveCyclicManyMany xs ys
:}
-}
convolveCyclicManyMany ::
   (Shape.C sh, Eq sh, Integral n, Class.Floating a) =>
   Array (sh, Shape.Cyclic n) a ->
   Array (sh, Shape.Cyclic n) a ->
   Array (sh, Shape.Cyclic n) a
convolveCyclicManyMany =
   getConvolve $
   Class.switchFloating
      (Convolve convolveCyclicManyManyReal)
      (Convolve convolveCyclicManyManyReal)
      (Convolve convolveCyclicManyManyComplex)
      (Convolve convolveCyclicManyManyComplex)


convolveCyclicReal ::
   (Integral n, Class.Real a) =>
   Array (Shape.Cyclic n) a ->
   Array (Shape.Cyclic n) a ->
   Array (Shape.Cyclic n) a
convolveCyclicReal x y =
   let Shape.Cyclic n = Array.shape x in
   let spectrum = Rank1.fourier Rank1.Forward $ Vector.zipComplex x y in
   let sqrHalf = (\v -> VectorSlice.mul v v) . takeShape (Spectrum.Half n) in
   Rank1.fourierCR $
   Vector.scale (0 :+ negate (recip $ 4 * fromIntegral n)) $
      sqrHalf spectrum
      |-|
      -- ToDo: optimize cyclicReverse with VectorSlice, we need only the first half
      -- or write it on a low-level
      -- or in the ST Monad
      Vector.conjugate (sqrHalf (Vector.cyclicReverse spectrum))

{- ToDo:
use joint two-real-as-one-complex transform
use Shape.NestedTuple Complex for reshaping array
Or does FFTW already perform this optimization?
-}
convolveCyclicOneManyReal ::
   (Shape.C sh, Integral n, Class.Real a) =>
   Array     (Shape.Cyclic n) a ->
   Array (sh, Shape.Cyclic n) a ->
   Array (sh, Shape.Cyclic n) a
convolveCyclicOneManyReal x y =
   Vector.scale (recip $ fromIntegral $ Shape.cyclicSize $ Array.shape x) $
   withExtraUnitDim Batch.fourierCR $
   Matrix.scaleColumns
      (Rank1.fourierRC x)
      (withExtraUnitDim Batch.fourierRC y)

convolveCyclicManyManyReal ::
   (Shape.C sh, Eq sh, Integral n, Class.Real a) =>
   Array (sh, Shape.Cyclic n) a ->
   Array (sh, Shape.Cyclic n) a ->
   Array (sh, Shape.Cyclic n) a
convolveCyclicManyManyReal x y =
   Vector.scale (recip $ fromIntegral $ Shape.cyclicSize $ Matrix.width x) $
   withExtraUnitDim Batch.fourierCR $
   Vector.mul
      (withExtraUnitDim Batch.fourierRC x)
      (withExtraUnitDim Batch.fourierRC y)


convolveCyclicComplex ::
   (Integral n, Class.Real a) =>
   Array (Shape.Cyclic n) (Complex a) ->
   Array (Shape.Cyclic n) (Complex a) ->
   Array (Shape.Cyclic n) (Complex a)
convolveCyclicComplex x y =
   Vector.scaleReal (recip $ fromIntegral $ Shape.cyclicSize $ Array.shape x) $
   Rank1.fourier Rank1.Backward $
   Vector.mul
      (Rank1.fourier Rank1.Forward x)
      (Rank1.fourier Rank1.Forward y)

convolveCyclicOneManyComplex ::
   (Shape.C sh, Integral n, Class.Real a) =>
   Array     (Shape.Cyclic n) (Complex a) ->
   Array (sh, Shape.Cyclic n) (Complex a) ->
   Array (sh, Shape.Cyclic n) (Complex a)
convolveCyclicOneManyComplex x y =
   Vector.scaleReal (recip $ fromIntegral $ Shape.cyclicSize $ Array.shape x) $
   Batch.fourier Rank1.Backward $
   Matrix.scaleColumns
      (Rank1.fourier Rank1.Forward x)
      (Batch.fourier Rank1.Forward y)

convolveCyclicManyManyComplex ::
   (Shape.C sh, Eq sh, Integral n, Class.Real a) =>
   Array (sh, Shape.Cyclic n) (Complex a) ->
   Array (sh, Shape.Cyclic n) (Complex a) ->
   Array (sh, Shape.Cyclic n) (Complex a)
convolveCyclicManyManyComplex x y =
   Vector.scaleReal
      (recip $ fromIntegral $ Shape.cyclicSize $ Matrix.width x) $
   Batch.fourier Rank1.Backward $
   Vector.mul
      (Batch.fourier Rank1.Forward x)
      (Batch.fourier Rank1.Forward y)



{- |
prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll (genVector number_) $ \y ->
   convolve x y =~= Naive.convolve x y
:}
-}
convolve ::
   (Integral n, Bits n, Class.Floating a) =>
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a
convolve =
   getConvolve $
   Class.switchFloating
      (Convolve convolveReal)
      (Convolve convolveReal)
      (Convolve convolveComplex)
      (Convolve convolveComplex)

{- ToDo:
make a variant returning VectorSlice
saves a copy for toCyclic in Rank1.convolve frontend.
-}
convolveReal ::
   (Integral n, Bits n, Class.Real a) =>
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a
convolveReal x y =
   let Shape.ZeroBased nx = Array.shape x in
   let Shape.ZeroBased ny = Array.shape y in
   if' (nx==0) x $ if' (ny==0) y $
   let n = max 1 (nx+ny) - 1 in
   -- ToDo: better use a ceilingEfficientEven
   let nsmooth = ceilingEfficient n in
   VectorSlice.toVector $ takeShape (Shape.ZeroBased n) $
   convolveCyclicReal (padCyclic nsmooth x) (padCyclic nsmooth y)


convolveComplex ::
   (Integral n, Bits n, Class.Real a) =>
   Array (Shape.ZeroBased n) (Complex a) ->
   Array (Shape.ZeroBased n) (Complex a) ->
   Array (Shape.ZeroBased n) (Complex a)
convolveComplex x y =
   let Shape.ZeroBased nx = Array.shape x in
   let Shape.ZeroBased ny = Array.shape y in
   if' (nx==0) x $ if' (ny==0) y $
   let n = max 1 (nx+ny) - 1 in
   let nsmooth = ceilingEfficient n in
   VectorSlice.scaleReal (recip $ fromIntegral nsmooth) $
   takeShape (Shape.ZeroBased n) $
   Rank1.fourier Rank1.Backward $
   -- ToDo: Fourier transforms could be done in a batch - would this be faster, at all?
   -- ToDo: do mul and scale in one go, however it's scaleReal
   Vector.mul
      (Rank1.fourier Rank1.Forward (padCyclic nsmooth x))
      (Rank1.fourier Rank1.Forward (padCyclic nsmooth y))


{- |
prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll genVectorShape $ \height ->
   QC.forAll genVectorShape $ \width ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \ys ->

   convolveOneMany x ys
   =~=
   Array2.fromRowArray
      (CShape.convolve (Array.shape x) width)
      (fmap (convolve x) (Array2.toRowArray ys))
:}

prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll genVectorShape $ \height ->
   QC.forAll genVectorShape $ \width ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \ys ->

   convolveOneMany x ys
   =~=
   Naive.convolveOneMany x ys
:}
-}
convolveOneMany ::
   (Shape.C sh, Eq sh, Integral n, Bits n, Class.Floating a) =>
   Array     (Shape.ZeroBased n) a ->
   Array (sh, Shape.ZeroBased n) a ->
   Array (sh, Shape.ZeroBased n) a
convolveOneMany x y =
   let shz@(Shape.ZeroBased n) =
         CShape.convolve (Array.shape x) (Matrix.width y) in
   let nsmooth = ceilingEfficient n in
   BlockMatrix.toMatrix $
   takeCenter (Shape.Zero ::+ shz ::+ Shape.ZeroBased (nsmooth-n)) $
   convolveOneManyPad nsmooth x y

{- |
Target size must be at least @max 0 (nx+ny-1)@.

prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll genVectorShape $ \height ->
   QC.forAll genVectorShape $ \width ->
   let Shape.ZeroBased n = CShape.convolve (Array.shape x) width in
   forChoose (n,n+100) $ \m ->
   QC.forAll (genVectorForShape number_ (height, width)) $ \ys ->

   convolveOneManyPad m x ys
   =~=
   padBlock (Shape.ZeroBased m) (Naive.convolveOneMany x ys)
:}
-}
convolveOneManyPad ::
   (Shape.C sh, Eq sh, Integral n, Class.Floating a) =>
   n ->
   Array     (Shape.ZeroBased n) a ->
   Array (sh, Shape.ZeroBased n) a ->
   Array (sh, Shape.ZeroBased n) a
convolveOneManyPad n x y =
   acyclicWidth $ convolveCyclicOneMany (padCyclic n x) (padBlockCyclic n y)

{- |
prop> :{
   QC.forAll genVectorShape $ \height ->
   QC.forAll genVectorShape $ \widthX ->
   QC.forAll genVectorShape $ \widthY ->
   QC.forAll (genVectorForShape number_ (height, widthX)) $ \xs ->
   QC.forAll (genVectorForShape number_ (height, widthY)) $ \ys ->

   let zs = convolveManyMany xs ys in
   zs
   =~=
   Array2.fromRowArray (Matrix.width zs)
      (BoxedArray.zipWith convolve
         (Array2.toRowArray xs) (Array2.toRowArray ys))
:}

prop> :{
   QC.forAll genVectorShape $ \height ->
   QC.forAll genVectorShape $ \widthX ->
   QC.forAll genVectorShape $ \widthY ->
   QC.forAll (genVectorForShape number_ (height, widthX)) $ \xs ->
   QC.forAll (genVectorForShape number_ (height, widthY)) $ \ys ->

   convolveManyMany xs ys
   =~=
   Naive.convolveManyMany xs ys
:}
-}
convolveManyMany ::
   (Shape.C sh, Eq sh, Integral n, Bits n, Class.Floating a) =>
   Array (sh, Shape.ZeroBased n) a ->
   Array (sh, Shape.ZeroBased n) a ->
   Array (sh, Shape.ZeroBased n) a
convolveManyMany x y =
   let shz@(Shape.ZeroBased n) =
         CShape.convolve (Matrix.width x) (Matrix.width y) in
   let nsmooth = ceilingEfficient n in
   BlockMatrix.toMatrix $
   takeCenter (Shape.Zero ::+ shz ::+ Shape.ZeroBased (nsmooth-n)) $
   convolveManyManyPad nsmooth x y

{- |
Target size must be at least @max 0 (nx+ny-1)@.

prop> :{
   QC.forAll genVectorShape $ \height ->
   QC.forAll genVectorShape $ \widthX ->
   QC.forAll genVectorShape $ \widthY ->
   QC.forAll (genVectorForShape number_ (height, widthX)) $ \xs ->
   QC.forAll (genVectorForShape number_ (height, widthY)) $ \ys ->
   let Shape.ZeroBased n = CShape.convolve widthX widthY in
   forChoose (n,n+100) $ \m ->

   convolveManyManyPad m xs ys
   =~=
   padBlock (Shape.ZeroBased m) (convolveManyMany xs ys)
:}
-}
convolveManyManyPad ::
   (Shape.C sh, Eq sh, Integral n, Class.Floating a) =>
   n ->
   Array (sh, Shape.ZeroBased n) a ->
   Array (sh, Shape.ZeroBased n) a ->
   Array (sh, Shape.ZeroBased n) a
convolveManyManyPad n x y =
   acyclicWidth $
   convolveCyclicManyMany (padBlockCyclic n x) (padBlockCyclic n y)



{-
ToDo:
rename to addOverlapping
-}
{-
ToDo:
openblas has sgeadd
mkl has mkl_somatadd
-}
{- |
ZeroBased size must be at most 2*stepSize.

prop> :{
   QC.forAll (genVector number_) $ \dat (QC.Positive width) ->
   forChoose (div (width+1) 2, width) $ \stepSize ->
   let (numBlocks,r) = divMod (Shape.zeroBasedSize $ Array.shape dat) width in
   let (xs,x) = Array.split $
         Array.reshape
            ((Shape.ZeroBased numBlocks,
               Shape.ZeroBased width) ::+ Shape.ZeroBased r) dat in
   forChoose (max r $ width-stepSize, width) $ \rr ->
   let completeSize = numBlocks*stepSize + rr in
   QC.counterexample ("numBlocks = " ++ show numBlocks) $
   QC.counterexample ("completeSize = " ++ show completeSize) $

   arrange completeSize stepSize xs x
   =~=
   arrangeViaMatrixAdd completeSize stepSize xs x
--   arrangeMutable completeSize stepSize xs x
:}
-}
arrange, arrangeMutable ::
   HasCallStack =>
--   (Show sh, Show n, Show a) =>
   (Shape.C sh, Integral n, Class.Floating a) =>
   n -> n ->
   Array (sh, Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a
arrange completeSize stepSize xs x =
--   trace ("arrange " ++ show (completeSize, stepSize, Array.shape xs, Array.shape x)) $
   let (height, width@(Shape.ZeroBased m)) = Array.shape xs in
   let n = Shape.size height in
   if' (n==0) (padShape (Shape.ZeroBased completeSize) x) $
   let height1 = Shape.ZeroBased $ pred n in
   let submat box =
         BlockMatrix.submatrix .
         Array.reshape
            (Subshape.focus (ViewMatrix.top . ViewMatrix.bottom) $
             Subshape.focus (ViewMatrix.left . ViewMatrix.right) $
             Subshape.submatrixFromMatrix box) in
   let z = Shape.Zero in
   let step = Shape.ZeroBased stepSize in
   let over = Shape.ZeroBased (m-stepSize) in
   let gap = Shape.ZeroBased (2*stepSize-m) in
   let a = submat (z ::+ height1 ::+ (), step ::+ over ::+ z) xs in
   let b = submat (() ::+ height1 ::+ z, z ::+ over ::+ step) xs in
   let c = submat (() ::+ height1 ::+ z, over ::+ gap ::+ over) xs in
--   trace ("a = " ++ show (BlockMatrix.toMatrix (BlockMatrix.submatrix a))) $
--   trace ("b = " ++ show (BlockMatrix.toMatrix (BlockMatrix.submatrix b))) $
--   trace ("c = " ++ show (BlockMatrix.toMatrix (BlockMatrix.submatrix c))) $
   let abc = BlockMatrix.toMatrix $
         BlockMatrix.toMatrix a |+| BlockMatrix.toMatrix b  &|||  c in
   let tailSize = completeSize - fromIntegral n * stepSize in
--   trace ("tailSize = " ++ show tailSize) $
   Array.reshape (Shape.ZeroBased completeSize) $ Chunk.toVector $
   Chunk.simple (takeShape (Shape.ZeroBased stepSize) xs)
   <>
   Chunk.simple (Array.mapShape (Shape.ZeroBased . Shape.size) abc)
   <>
   Chunk.simple
      (padShape (Shape.ZeroBased tailSize)
         (VectorSlice.sliceVector (Slice.sub stepSize tailSize . Slice.right)
            (Array.reshape ((height1,width) ::+ width) xs))
       |+|
       padShape (Shape.ZeroBased tailSize) x)

arrangeMutable completeSize stepSize xs x = MutArray.runST (do
   let addAt arr k row =
         MutVector.add
            (MutVector.sliceVector
               (Slice.sub k (Shape.zeroBasedSize $ VectorSlice.shape row)) arr)
            row
   arr <- MutVector.new (Shape.ZeroBased completeSize) Scalar.zero
   zipWithM_ (addAt arr)
      (iterate (stepSize+) 0)
      (VectorSlice.slicesVector Slice.rows xs)
   addAt arr (stepSize * fromIntegral (Shape.size (Matrix.height xs))) x
   return arr)

arrangeViaMatrixAdd ::
--   HasCallStack =>
--   (Show sh, Show n, Show a) =>
   (Shape.C sh, Eq sh, Integral n, Class.Floating a) =>
   n -> n ->
   Array (sh, Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a
-- ToDo: initialize with copied submatrix and only fill remainder with zero
arrangeViaMatrixAdd completeSize stepSize xs x = MutArray.runST (do
--   ys <- MutVector.new (Shape.ZeroBased completeSize) Scalar.zero
   ys <- MutArray.new (Shape.ZeroBased completeSize) Scalar.zero
   let xssplit =
--         Array.mapShape Subshape.submatrixFromMatrix $
         Matrix.mapWidth (splitShape (Shape.ZeroBased stepSize)) xs
   let subshape =
         Subshape.Cons {
            Subshape.start = 0,
            Subshape.size = fromIntegral completeSize,
            Subshape.layout = Layout.Submatrix $ fromIntegral stepSize,
            Subshape.shape = (Matrix.height xs, Shape.ZeroBased stepSize)
         }
   BlockMatrix.toMutableSubmatrix
      (BlockMatrix.submatrix $
         Array.mapShape
            (Subshape.focus ViewMatrix.left . Subshape.submatrixFromMatrix)
            xssplit)
      (MutArray.reshape subshape ys)
{-
      (MutArray.mapShape
         ((\shape -> shape{Subshape.shape = (Matrix.height y, Array.shape x)}) .
            ViewMatrix.submatrixFromMatrix)
         ys)
-}
   let (_ ::+ overlapSize) = Matrix.width xssplit
   MutMatrix.addSubmatrix
      Scalar.one
      (Matrix.viewMatrix ViewMatrix.right xssplit)
      Scalar.one
      (MutMatrix.wrap $
       -- move window one row down
       MutArray.reshape
         (subshape {
            Subshape.start = fromIntegral stepSize,
            Subshape.shape = (Matrix.height xs, overlapSize)
         })
         ys)
   let addAt k row =
         MutVector.add
            (MutVector.sliceVector
               (Slice.sub k (Shape.zeroBasedSize $ VectorSlice.shape row)) ys)
            row
   addAt (stepSize * fromIntegral (Shape.size (Matrix.height xs))) x
   return ys)


{- |
prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll (genVector number_) $ \y ->
   convolve x y =~= convolveShortLong x y
:}

prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll (genVector number_) $ \y0 ->
   QC.forAll (genVector number_) $ \y1 ->
   QC.forAll (genVector number_) $ \y2 ->
   let y = y0<>y1<>y2 in
   convolve x y =~= convolveShortLong x y
:}
-}
convolveShortLong ::
--   (Show n, Show a) =>
   (Integral n, Bits n, Class.Floating a) =>
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a
convolveShortLong x y =
   let Shape.ZeroBased nx = Array.shape x in
   let Shape.ZeroBased ny = Array.shape y in
   let n = max 1 (nx+ny) - 1 in
   let nsmooth = ceilingPowerOfTwo (max 1 (2*nx)) in
   let chunkSize = nsmooth+1-nx in
   let (numCompleteChunks, tlSize) = divMod ny chunkSize in
   let (block,tl) =
         Array.split $
         Array.reshape
            ((Shape.ZeroBased numCompleteChunks, Shape.ZeroBased chunkSize)
               ::+ Shape.ZeroBased tlSize)
            y in
   {-
   We call convolveOneManyPad not for extra padding,
   but for avoiding takeCenter.
   -}
   arrangeViaMatrixAdd n chunkSize (convolveOneManyPad nsmooth x block) (convolve x tl)

convolveShortLongInterval ::
   (Show n, Show a) =>
   (Integral n, Bits n, Class.Floating a) =>
   Array (Shape.Interval n) a ->
   Array (Shape.Interval n) a ->
   Array (Shape.Interval n) a
convolveShortLongInterval x y =
   let (Shape.Interval xFrom xTo) = Array.shape x in
   let (Shape.Interval yFrom yTo) = Array.shape y in
   Array.reshape (Shape.Interval (xFrom+yFrom) (xTo+yTo)) $
   convolveShortLong
      (Array.reshape (Shape.ZeroBased (xTo-xFrom+1)) x)
      (Array.reshape (Shape.ZeroBased (yTo-yFrom+1)) y)

{-
ToDo:
The implementation can be optimized
such that outer elements are not computed.
They are dropped anyway.
Similar to convolveLongToShort.
However, how much can we save with this optimization?
The number of small convolutions remains the same.
We could at least save something at the addOverlap stage.
-}
convolveShortLongCut ::
   (Show n, Show a) =>
   (Integral n, Bits n, Class.Floating a) =>
   Array (Shape.Interval n) a ->
   Array (Shape.Interval n) a ->
   VectorSlice.T Layout.Slice (Shape.Interval n) a
convolveShortLongCut x y =
   let (Shape.Interval xFrom xTo) = Array.shape x in
   let (Shape.Interval yFrom yTo) = Array.shape y in
   VectorSlice.sliceVector
      (Slice.subInterval (Shape.Interval (yFrom+xTo) (xFrom+yTo))) $
   convolveShortLongInterval x y


-- ToDo: the shortening behavior better fits a correlate function
{- |
Convolve two long vectors and keep only the short time of maximal overlap.
For a non-empty result, @y@ must be at least as long @x@.

This requires only computation time proportional to
of @(size x + size y) * log (size y - size x)@.

prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll (genVector number_) $ \y ->

   convolveLongToShort x y
   =~=
   Naive.convolveLongToShort x y
:}

prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll (genVector number_) $ \y ->

   convolveLongToShort x y
   =~=
   (Vector.drop (pred $ Shape.zeroBasedSize $ Array.shape x) $
    Vector.take (Shape.zeroBasedSize $ Array.shape y) $
      convolve x y)
:}

prop> :{
   QC.forAll (genVector number_) $ \x ->
   QC.forAll (genVector number_) $ \y ->
   forChoose (0, Shape.zeroBasedSize (Array.shape y)) $ \m ->

   convolveLongToShort x (Vector.drop m y)
   =~=
   Vector.drop m (convolveLongToShort x y)
:}

prop> :{
   QC.forAll (genVector number_) $ \x ->
   let n = Shape.zeroBasedSize (Array.shape x) in
   QC.forAll (genVectorForShape number_ (Array.shape x)) $ \y0 ->
   QC.forAll (genVector number_) $ \y1 ->
   forChoose (0, Shape.zeroBasedSize (Array.shape y1)) $ \m ->
   let y = y0<>y1 in

   convolveLongToShort x (Vector.take (n+m) y)
   =~=
   Vector.take (m+1) (convolveLongToShort x y)
:}
-}
{-
   ABCDEF * 12345678 =
        12345678
       12345678
      12345678
     12345678
    12345678
   12345678
-}
convolveLongToShort ::
--   (Show n, Show a) =>
   (Integral n, Bits n, Class.Floating a) =>
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a ->
   Array (Shape.ZeroBased n) a
convolveLongToShort x y =
   let Shape.ZeroBased nx = Array.shape x in
   let Shape.ZeroBased ny = Array.shape y in
   let overlap = max 0 $ 1+ny-nx in
   if' (overlap == 0) (Vector.zero (Shape.ZeroBased 0)) $
   if' (overlap == 1)
         (Vector.constant (Shape.ZeroBased 1) $
          Vector.dot x $ Vector.reverse y) $
   let nsmooth = ceilingPowerOfTwo (2*overlap) in
   let stepSize = nsmooth - overlap + 1 in
   let numChunks = div nx stepSize in
   let x0size = nx - numChunks * stepSize in
   -- ToDo: save the copy into x0 and xBlock with VectorSlice and BlockMatrix
   let (x0,xBlock) =
         Array.split $
         Array.reshape
            (Shape.ZeroBased x0size
               ::+ (Shape.ZeroBased numChunks, Shape.ZeroBased stepSize)) x in
   let yRows =
         fmap (\k ->
            VectorSlice.sliceVector
               (Slice.reshape (Shape.Cyclic nsmooth) .
                Slice.sub ((numChunks-1-k)*stepSize) nsmooth)
               y) $
         BoxedArray.indices (Shape.ZeroBased numChunks) in
   let y0 = VectorSlice.sliceVector (Slice.drop (numChunks*stepSize)) y in
   (if x0size==0
      then id
      else
         let z0size =
               ceilingEfficient $ x0size+Shape.zeroBasedSize(Array.shape y)-1 in
         VectorSlice.add
            (VectorSlice.sliceVector (Slice.sub (x0size-1) overlap) $ acyclic $
             convolveCyclic (padCyclic z0size x0) (padCyclic z0size y0))) $
   Matrix.sumColumns $
   Matrix.takeRight $
   Matrix.reshapeWidth
      (Shape.ZeroBased (nsmooth-overlap) ::+ Shape.ZeroBased overlap) $
   convolveCyclicManyMany
      (padBlockCyclic nsmooth xBlock)
      (Matrix.fromRowArray (Shape.Cyclic nsmooth) yRows)

{- |
Attention:
It does not hold:

fromZeroBased convolveLongToShort x y
==
convolveLongToShortInterval (fromZeroBased x) (fromZeroBased y)

because convolveLongToShortInterval returns a shape with a non-zero offset.

The identity will only be true for a correlation function.
-}
convolveLongToShortInterval ::
--   (Show n, Show a) =>
   (Integral n, Bits n, Class.Floating a) =>
   Array (Shape.Interval n) a ->
   Array (Shape.Interval n) a ->
   Array (Shape.Interval n) a
convolveLongToShortInterval x y =
   Array.reshape
      (CShape.correlate
         (Shape.intervalNegate $ Array.shape x) (Array.shape y)) $
   convolveLongToShort
      (Array.mapShape CShape.unshift x)
      (Array.mapShape CShape.unshift y)