content: {-# LANGUAGE TypeFamilies #-}
  {-# LANGUAGE FlexibleContexts #-}
  {-# LANGUAGE FlexibleInstances #-}
  {-# LANGUAGE MultiParamTypeClasses #-}
  {-# LANGUAGE FunctionalDependencies #-}
  {-# LANGUAGE UndecidableInstances #-}
  
  -----------------------------------------------------------------------------
  -- |
  -- Module      :  Numeric.Container
  -- Copyright   :  (c) Alberto Ruiz 2010
  -- License     :  GPL-style
  --
  -- Maintainer  :  Alberto Ruiz <aruiz@um.es>
  -- Stability   :  provisional
  -- Portability :  portable
  --
  -- 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 Numeric.Container (
      -- * Basic functions
      module Data.Packed,
      constant, linspace,
      diag, ident,
      ctrans,
      -- * Generic operations
      Container(..),
      -- * Matrix product
      Product(..),
      optimiseMult,
      mXm,mXv,vXm,(<.>),Mul(..),LSDiv(..),
      outer, kronecker,
      -- * Random numbers
      RandDist(..),
      randomVector,
      gaussianSample,
      uniformSample,
      meanCov,
      -- * Element conversion
      Convert(..),
      Complexable(),
      RealElement(),
  
      RealOf, ComplexOf, SingleOf, DoubleOf,
  
      IndexOf,
      module Data.Complex,
      -- * Input / Output
      dispf, disps, dispcf, vecdisp, latexFormat, format,
      loadMatrix, saveMatrix, fromFile, fileDimensions,
      readMatrix,
      fscanfVector, fprintfVector, freadVector, fwriteVector,
      -- * Experimental
      build', konst'
  ) where
  
  import Data.Packed
  import Data.Packed.Internal(constantD)
  import Numeric.ContainerBoot
  import Numeric.Chain
  import Numeric.IO
  import Data.Complex
  import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
  import Data.Packed.Random
  
  ------------------------------------------------------------------
  
  {- | creates a vector with a given number of equal components:
  
  @> constant 2 7
  7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
  -}
  constant :: Element a => a -> Int -> Vector a
  -- constant x n = runSTVector (newVector x n)
  constant = constantD-- about 2x faster
  
  {- | Creates a real vector containing a range of values:
  
  @\> linspace 5 (-3,7)
  5 |> [-3.0,-0.5,2.0,4.5,7.0]@
  
  Logarithmic spacing can be defined as follows:
  
  @logspace n (a,b) = 10 ** linspace n (a,b)@
  -}
  linspace :: (Enum e, Container Vector e) => Int -> (e, e) -> Vector e
  linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1]
      where s = (b-a)/fromIntegral (n-1)
  
  -- | Dot product: @u \<.\> v = dot u v@
  (<.>) :: Product t => Vector t -> Vector t -> t
  infixl 7 <.>
  (<.>) = dot
  
  
  
  --------------------------------------------------------
  
  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
      (<>) = mXm
  
  instance Mul Matrix Vector Vector where
      (<>) m v = flatten $ m <> asColumn v
  
  instance Mul Vector Matrix Vector where
      (<>) v m = flatten $ asRow v <> m
  
  --------------------------------------------------------
  
  class LSDiv b c | b -> c, c->b where
   infixl 7 <\>
   -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
   (<\>)  :: Field t => Matrix t -> b t -> c t
  
  instance LSDiv Vector Vector where
      m <\> v = flatten (linearSolveSVD m (reshape 1 v))
  
  instance LSDiv Matrix Matrix where
      (<\>) = linearSolveSVD
  
  --------------------------------------------------------
  
  -- | Compute mean vector and covariance matrix of the rows of a matrix.
  meanCov :: Matrix Double -> (Vector Double, Matrix Double)
  meanCov x = (med,cov) where
      r    = rows x
      k    = 1 / fromIntegral r
      med  = konst k r `vXm` x
      meds = konst 1 r `outer` med
      xc   = x `sub` meds
      cov  = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)
  
  