    {-# LANGUAGE CPP #-}
    {-# LANGUAGE FlexibleContexts #-}
    {-# LANGUAGE FlexibleInstances #-}
    {-# LANGUAGE MultiParamTypeClasses #-}
    {-# LANGUAGE FunctionalDependencies #-}

    -----------------------------------------------------------------------------
    -- |
    -- Module      :  Internal.Container
    -- Copyright   :  (c) Alberto Ruiz 2010-14
    -- License     :  BSD3
    -- Maintainer  :  Alberto Ruiz
    -- Stability   :  provisional
    --
    -- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
    --
    -- The 'Container' class is used to define optimized generic functions which work
    -- on 'Vector' and 'Matrix' with real or complex elements.
    --
    -- Some of these functions are also available in the instances of the standard
    -- numeric Haskell classes provided by "Numeric.LinearAlgebra".
    --
    -----------------------------------------------------------------------------

    module Internal.Container where

    import Internal.Vector
    import Internal.Matrix
    import Internal.Element
    import Internal.Numeric
    import Internal.Algorithms(Field,linearSolveSVD,Herm,mTm)
    #if MIN_VERSION_base(4,11,0)
    import Prelude hiding ((<>))
    #endif
    ------------------------------------------------------------------

    {- | Creates a real vector containing a range of values:

    >>> linspace 5 (-3,7::Double)
    [-3.0,-0.5,2.0,4.5,7.0]
    it :: Vector Double

    >>> linspace 5 (8,3:+2) :: Vector (Complex Double)
    [8.0 :+ 0.0,6.75 :+ 0.5,5.5 :+ 1.0,4.25 :+ 1.5,3.0 :+ 2.0]
    it :: Vector (Complex Double)

    Logarithmic spacing can be defined as follows:

    @logspace n (a,b) = 10 ** linspace n (a,b)@
    -}
    linspace :: (Fractional e, Container Vector e) => Int -> (e, e) -> Vector e
    linspace :: forall e.
    (Fractional e, Container Vector e) =>
    Int -> (e, e) -> Vector e
    linspace Int
    0 (e, e)
    _     = [e] -> Vector e
    forall a. Storable a => [a] -> Vector a
    fromList[]
    linspace Int
    1 (e
    a,e
    b) = [e] -> Vector e
    forall a. Storable a => [a] -> Vector a
    fromList[(e
    ae -> e -> e
    forall a. Num a => a -> a -> a
    +e
    b)e -> e -> e
    forall a. Fractional a => a -> a -> a
    /e
    2]
    linspace Int
    n (e
    a,e
    b) = e -> Vector e -> Vector e
    forall (c :: * -> *) e. Container c e => e -> c e -> c e
    addConstant e
    a (Vector e -> Vector e) -> Vector e -> Vector e
    forall a b. (a -> b) -> a -> b
    $ e -> Vector e -> Vector e
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale e
    s (Vector e -> Vector e) -> Vector e -> Vector e
    forall a b. (a -> b) -> a -> b
    $ [e] -> Vector e
    forall a. Storable a => [a] -> Vector a
    fromList ([e] -> Vector e) -> [e] -> Vector e
    forall a b. (a -> b) -> a -> b
    $ (Int -> e) -> [Int] -> [e]
    forall a b. (a -> b) -> [a] -> [b]
    map Int -> e
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral [Int
    0 .. Int
    nInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1]
        where s :: e
    s = (e
    be -> e -> e
    forall a. Num a => a -> a -> a
    -e
    a)e -> e -> e
    forall a. Fractional a => a -> a -> a
    /Int -> e
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Int
    nInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1)

    --------------------------------------------------------------------------------

    infixr 8 <.>
    {- | An infix synonym for 'dot'

    >>> vector [1,2,3,4] <.> vector [-2,0,1,1]
    5.0

    >>> let 𝑖 = 0:+1 :: C
    >>> fromList [1+𝑖,1] <.> fromList [1,1+𝑖]
    2.0 :+ 0.0

    -}

    (<.>) :: Numeric t => Vector t -> Vector t -> t
    <.> :: forall t. Numeric t => Vector t -> Vector t -> t
    (<.>) = Vector t -> Vector t -> t
    forall t. Numeric t => Vector t -> Vector t -> t
    dot

    {- | dense matrix-vector product

    >>> let m = (2><3) [1..]
    >>> m
    (2><3)
     [ 1.0, 2.0, 3.0
     , 4.0, 5.0, 6.0 ]

    >>> let v = vector [10,20,30]

    >>> m #> v
    [140.0,320.0]
    it :: Vector Numeric.LinearAlgebra.Data.R

    -}
    infixr 8 #>
    (#>) :: Numeric t => Matrix t -> Vector t -> Vector t
    #> :: forall t. Numeric t => Matrix t -> Vector t -> Vector t
    (#>) = Matrix t -> Vector t -> Vector t
    forall t. Product t => Matrix t -> Vector t -> Vector t
    mXv

    -- | dense matrix-vector product
    app :: Numeric t => Matrix t -> Vector t -> Vector t
    app :: forall t. Numeric t => Matrix t -> Vector t -> Vector t
    app = Matrix t -> Vector t -> Vector t
    forall t. Numeric t => Matrix t -> Vector t -> Vector t
    (#>)

    infixl 8 <#
    -- | dense vector-matrix product
    (<#) :: Numeric t => Vector t -> Matrix t -> Vector t
    <# :: forall t. Numeric t => Vector t -> Matrix t -> Vector t
    (<#) = Vector t -> Matrix t -> Vector t
    forall t. Product t => Vector t -> Matrix t -> Vector t
    vXm

    --------------------------------------------------------------------------------

    class Mul a b c | a b -> c where
     infixl 7 <>
     -- | Matrix-matrix, matrix-vector, and vector-matrix products.
     (<>)  :: Product t => a t -> b t -> c t

    instance Mul Matrix Matrix Matrix where
        <> :: forall t. Product t => Matrix t -> Matrix t -> Matrix t
    (<>) = Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    mXm

    instance Mul Matrix Vector Vector where
        <> :: forall t. Product t => Matrix t -> Vector t -> Vector t
    (<>) Matrix t
    m Vector t
    v = Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten (Matrix t -> Vector t) -> Matrix t -> Vector t
    forall a b. (a -> b) -> a -> b
    $ Matrix t
    m Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Vector t -> Matrix t
    forall a. Storable a => Vector a -> Matrix a
    asColumn Vector t
    v

    instance Mul Vector Matrix Vector where
        <> :: forall t. Product t => Vector t -> Matrix t -> Vector t
    (<>) Vector t
    v Matrix t
    m = Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten (Matrix t -> Vector t) -> Matrix t -> Vector t
    forall a b. (a -> b) -> a -> b
    $ Vector t -> Matrix t
    forall a. Storable a => Vector a -> Matrix a
    asRow Vector t
    v Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix t
    m

    --------------------------------------------------------------------------------

    {- | Least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)

    @
    a = (3><2)
     [ 1.0,  2.0
     , 2.0,  4.0
     , 2.0, -1.0 ]
    @

    @
    v = vector [13.0,27.0,1.0]
    @

    >>> let x = a <\> v
    >>> x
    [3.0799999999999996,5.159999999999999]
    it :: Vector Numeric.LinearAlgebra.Data.R

    >>> a #> x
    [13.399999999999999,26.799999999999997,0.9999999999999991]
    it :: Vector Numeric.LinearAlgebra.Data.R

    It also admits multiple right-hand sides stored as columns in a matrix.

    -}
    infixl 7 <\>
    (<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t
    <\> :: forall (c :: * -> *) t.
    (LSDiv c, Field t) =>
    Matrix t -> c t -> c t
    (<\>) = Matrix t -> c t -> c t
    forall t. Field t => Matrix t -> c t -> c t
    forall (c :: * -> *) t.
    (LSDiv c, Field t) =>
    Matrix t -> c t -> c t
    linSolve

    class LSDiv c
      where
        linSolve :: Field t => Matrix t -> c t -> c t

    instance LSDiv Vector
      where
        linSolve :: forall t. Field t => Matrix t -> Vector t -> Vector t
    linSolve Matrix t
    m Vector t
    v = Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten (Matrix t -> Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolveSVD Matrix t
    m (Int -> Vector t -> Matrix t
    forall t. Storable t => Int -> Vector t -> Matrix t
    reshape Int
    1 Vector t
    v))

    instance LSDiv Matrix
      where
        linSolve :: forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linSolve = Matrix t -> Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolveSVD

    --------------------------------------------------------------------------------

    class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
      where
        -- |
        -- >>> build 5 (**2) :: Vector Double
        -- [0.0,1.0,4.0,9.0,16.0]
        -- it :: Vector Double
        --
        -- Hilbert matrix of order N:
        --
        -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
        -- >>> putStr . dispf 2 $ hilb 3
        -- 3x3
        -- 1.00  0.50  0.33
        -- 0.50  0.33  0.25
        -- 0.33  0.25  0.20
        --
        build :: d -> f -> c e

    instance Container Vector e => Build Int (e -> e) Vector e
      where
        build :: Int -> (e -> e) -> Vector e
    build = Int -> (e -> e) -> Vector e
    IndexOf Vector -> ArgOf Vector e -> Vector e
    forall (c :: * -> *) e.
    Container c e =>
    IndexOf c -> ArgOf c e -> c e
    build'

    instance (Num e, Container Vector e) => Build (Int,Int) (e -> e -> e) Matrix e
      where
        build :: (Int, Int) -> (e -> e -> e) -> Matrix e
    build = (Int, Int) -> (e -> e -> e) -> Matrix e
    IndexOf Matrix -> ArgOf Matrix e -> Matrix e
    forall (c :: * -> *) e.
    Container c e =>
    IndexOf c -> ArgOf c e -> c e
    build'

    --------------------------------------------------------------------------------

    -- @dot u v = 'udot' ('conj' u) v@
    dot :: (Numeric t) => Vector t -> Vector t -> t
    dot :: forall t. Numeric t => Vector t -> Vector t -> t
    dot Vector t
    u Vector t
    v = Vector t -> Vector t -> t
    forall e. Product e => Vector e -> Vector e -> e
    udot (Vector t -> Vector t
    forall (c :: * -> *) e. Container c e => c e -> c e
    conj Vector t
    u) Vector t
    v

    --------------------------------------------------------------------------------

    optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
    optimiseMult :: forall t. Monoid (Matrix t) => [Matrix t] -> Matrix t
    optimiseMult = [Matrix t] -> Matrix t
    forall a. Monoid a => [a] -> a
    mconcat

    --------------------------------------------------------------------------------

    {- | Compute mean vector and covariance matrix of the rows of a matrix.

    >>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (trustSym $ diagl [2,3])
    ([3.9933155655086696,5.061409102770331],Herm (2><2)
     [    1.9963242906624408, -4.227815571404954e-2
     , -4.227815571404954e-2,    3.2003833097832857 ])
    it :: (Vector Double, Herm Double)
    -}
    meanCov :: Matrix Double -> (Vector Double, Herm Double)
    meanCov :: Matrix Double -> (Vector Double, Herm Double)
    meanCov Matrix Double
    x = (Vector Double
    med,Herm Double
    cov) where
        r :: Int
    r    = Matrix Double -> Int
    forall t. Matrix t -> Int
    rows Matrix Double
    x
        k :: Double
    k    = Double
    1 Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Int
    r
        med :: Vector Double
    med  = Double -> Int -> Vector Double
    forall e d (c :: * -> *). Konst e d c => e -> d -> c e
    konst Double
    k Int
    r Vector Double -> Matrix Double -> Vector Double
    forall t. Product t => Vector t -> Matrix t -> Vector t
    `vXm` Matrix Double
    x
        meds :: Matrix Double
    meds = Double -> Int -> Vector Double
    forall e d (c :: * -> *). Konst e d c => e -> d -> c e
    konst Double
    1 Int
    r Vector Double -> Vector Double -> Matrix Double
    forall t. Product t => Vector t -> Vector t -> Matrix t
    `outer` Vector Double
    med
        xc :: Matrix Double
    xc   = Matrix Double
    x Matrix Double -> Matrix Double -> Matrix Double
    forall (c :: * -> *) e. Container c e => c e -> c e -> c e
    `sub` Matrix Double
    meds
        cov :: Herm Double
    cov  = Double -> Herm Double -> Herm Double
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale (Double -> Double
    forall a. Fractional a => a -> a
    recip (Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1))) (Matrix Double -> Herm Double
    forall t. Numeric t => Matrix t -> Herm t
    mTm Matrix Double
    xc)

    --------------------------------------------------------------------------------

    sortVector :: (Ord t, Element t) => Vector t -> Vector t
    sortVector :: forall t. (Ord t, Element t) => Vector t -> Vector t
    sortVector = Vector t -> Vector t
    forall a. (Element a, Ord a) => Vector a -> Vector a
    sortV

    {- |

    >>> m <- randn 4 10
    >>> disp 2 m
    4x10
    -0.31   0.41   0.43  -0.19  -0.17  -0.23  -0.17  -1.04  -0.07  -1.24
     0.26   0.19   0.14   0.83  -1.54  -0.09   0.37  -0.63   0.71  -0.50
    -0.11  -0.10  -1.29  -1.40  -1.04  -0.89  -0.68   0.35  -1.46   1.86
     1.04  -0.29   0.19  -0.75  -2.20  -0.01   1.06   0.11  -2.09  -1.58

    >>> disp 2 $ m ?? (All, Pos $ sortIndex (m!1))
    4x10
    -0.17  -1.04  -1.24  -0.23   0.43   0.41  -0.31  -0.17  -0.07  -0.19
    -1.54  -0.63  -0.50  -0.09   0.14   0.19   0.26   0.37   0.71   0.83
    -1.04   0.35   1.86  -0.89  -1.29  -0.10  -0.11  -0.68  -1.46  -1.40
    -2.20   0.11  -1.58  -0.01   0.19  -0.29   1.04   1.06  -2.09  -0.75

    -}
    sortIndex :: (Ord t, Element t) => Vector t -> Vector I
    sortIndex :: forall t. (Ord t, Element t) => Vector t -> Vector I
    sortIndex = Vector t -> Vector I
    forall a. (Element a, Ord a) => Vector a -> Vector I
    sortI

    ccompare :: (Ord t, Container c t) => c t -> c t -> c I
    ccompare :: forall t (c :: * -> *). (Ord t, Container c t) => c t -> c t -> c I
    ccompare = c t -> c t -> c I
    forall (c :: * -> *) e. (Container c e, Ord e) => c e -> c e -> c I
    ccompare'

    cselect :: (Container c t) => c I -> c t -> c t -> c t -> c t
    cselect :: forall (c :: * -> *) t.
    Container c t =>
    c I -> c t -> c t -> c t -> c t
    cselect = c I -> c t -> c t -> c t -> c t
    forall (c :: * -> *) t.
    Container c t =>
    c I -> c t -> c t -> c t -> c t
    cselect'

    {- | Extract elements from positions given in matrices of rows and columns.

    >>> r
    (3><3)
     [ 1, 1, 1
     , 1, 2, 2
     , 1, 2, 3 ]
    >>> c
    (3><3)
     [ 0, 1, 5
     , 2, 2, 1
     , 4, 4, 1 ]
    >>> m
    (4><6)
     [  0,  1,  2,  3,  4,  5
     ,  6,  7,  8,  9, 10, 11
     , 12, 13, 14, 15, 16, 17
     , 18, 19, 20, 21, 22, 23 ]

    >>> remap r c m
    (3><3)
     [  6,  7, 11
     ,  8, 14, 13
     , 10, 16, 19 ]

    The indexes are autoconformable.

    >>> c'
    (3><1)
     [ 1
     , 2
     , 4 ]
    >>> remap r c' m
    (3><3)
     [  7,  7,  7
     ,  8, 14, 14
     , 10, 16, 22 ]

    -}
    remap :: Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
    remap :: forall t. Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
    remap Matrix I
    i Matrix I
    j Matrix t
    m
        | Matrix I -> I
    forall (c :: * -> *) e. Container c e => c e -> e
    minElement Matrix I
    i I -> I -> Bool
    forall a. Ord a => a -> a -> Bool
    >= I
    0 Bool -> Bool -> Bool
    && Matrix I -> I
    forall (c :: * -> *) e. Container c e => c e -> e
    maxElement Matrix I
    i I -> I -> Bool
    forall a. Ord a => a -> a -> Bool
    < Int -> I
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m) Bool -> Bool -> Bool
    &&
          Matrix I -> I
    forall (c :: * -> *) e. Container c e => c e -> e
    minElement Matrix I
    j I -> I -> Bool
    forall a. Ord a => a -> a -> Bool
    >= I
    0 Bool -> Bool -> Bool
    && Matrix I -> I
    forall (c :: * -> *) e. Container c e => c e -> e
    maxElement Matrix I
    j I -> I -> Bool
    forall a. Ord a => a -> a -> Bool
    < Int -> I
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m) = Matrix I -> Matrix I -> Matrix t -> Matrix t
    forall t. Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
    remapM Matrix I
    i' Matrix I
    j' Matrix t
    m
        | Bool
    otherwise = [Char] -> Matrix t
    forall a. HasCallStack => [Char] -> a
    error ([Char] -> Matrix t) -> [Char] -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ [Char]
    "out of range index in remap"
      where
        [Matrix I
    i',Matrix I
    j'] = [Matrix I] -> [Matrix I]
    forall t. Element t => [Matrix t] -> [Matrix t]
    conformMs [Matrix I
    i,Matrix I
    j]