-- While working on this module you are encouraged to remove it and fix
-- any warnings in the module. See
--     http://hackage.haskell.org/trac/ghc/wiki/WorkingConventions#Warnings
-- for details  

-----------------------------------------------------------------------------
-- |
-- Module      :  Matrix
-- Copyright   :  (c) Lennart Schmitt
-- License     :  BSD-style (see the file libraries/base/LICENSE)
-- 
-- Maintainer  :  lennart...schmitt@<nospam>gmail.com
-- Stability   :  experimental
-- Portability :  portable
--
-- This module includes a few (standard) functions to work with matrixes.
-----------------------------------------------------------------------------
module Numeric.Matrix (
             -- * Data Types
             Matrix, -- | Represents a matrix
             RawMatrix, 
             Element, -- | The elements of a matrix

             -- * Convert a Matrix into another Type
             flatten,   -- | Flattens a matrix to a vector
             toColumns, -- | Converts the representation of a matrix to a list of vectors by its columns.
             toLists,   -- | Converts the representation of a matrix to a list of lists. (NOT equivalent with toRows!!!)
             toRows,    -- | Converts the representation of a matrix to a list of vectors by its rows.

             -- * Convert some Data into Matrix-Type
             asRow,     -- | Converts the representation of a matrix from a list of vectors to a matrix (by rows)
             fromLists, -- | Converts the representation of a matrix from a list of lists to a matrix (NOT equivalent with asRows!!!)
             fromListToQuadraticMatrix,

             -- * Matrix calculations/transformations
             (@@>), -- | Get the matrix-element in row x and col y (like "(!!)" for lists)
             eigenvalue,
             eigenvector,
             inv, -- | Calculates the invariant matrix of the input matrix
             identityMatrix,
             mapMatrix,
             reduceMatrix,
             scalarMultiplication,
             subtractMatrix,
             trans, -- | Transposes a matrix
             zipAllWith
       ) where

-- This module "extends" an existing matrix-module by using its datatypes and
-- it also uses linear algebra algorithms (so it is based on that modules too)
import Data.Packed.Matrix (Matrix, buildMatrix, fromLists, toLists, flatten, reshape, trans, rows,asRow,toRows, toColumns, (@@>), Element)
import Numeric.LinearAlgebra.Algorithms (eig, eigenvalues,inv)
-- There are also a few additional modules that are needed to work with the matrixes
import Numeric.Vector (Vector,RawVector,toList,fromList,count,maximumBy,transpose)
import Foreign.Marshal.Utils (fromBool)
import Data.Complex (realPart)

-- | A matrix represented by a list of lists.
type RawMatrix a = [[a]]

-- | Calculates the identity matrix (n x n) by given scale (n)
identityMatrix :: Int -> Matrix Double
identityMatrix i = buildMatrix i i (\ (x,y) -> fromBool (x == y))

-- | A simple map-Function which maps a given function on every element of the given matrix
mapMatrix :: (Double -> Double) -> Matrix Double -> Matrix Double
mapMatrix f m = fromLists.(map$map f).toLists $ m

-- | Calculates the scalarproduct (with a scalar and matrix)
scalarMultiplication :: Double -> Matrix Double -> Matrix Double
scalarMultiplication x m = mapMatrix ((*) x) m

-- | Calculates the difference (matrix) between two matrixes
subtractMatrix :: Matrix Double -> Matrix Double -> Matrix Double
subtractMatrix a b = fromListToQuadraticMatrix $ zipWith (-) (toList . flatten $a) (toList . flatten $b)

-- | Builds a quadratic matrix out of a list
fromListToQuadraticMatrix :: [Double] -> Matrix Double
fromListToQuadraticMatrix xs = (reshape (round$ sqrt(count xs)) (fromList xs))

-- | Calculates the eigenvalue of a matrix, e.g.
-- 
-- > eigenvalue (fromLists [[0.77143,-0.25714],[-0.42245,0.14082]]) 
-- 
-- returns 
-- 
-- > 0.9122456375784809 
eigenvalue :: Matrix Double -> Double
eigenvalue = ((maximumBy (compare . abs)) . (map realPart) . toList . eigenvalues)

-- | Calculates one eigenvector of a given matrix, e.g.
-- 
-- > eigenvector (fromLists [[-0.14081563757848092,-0.25714],[-0.42245,-0.7714256375784809]])
--
-- returns
--
-- > [0.8770950095147589,-0.48031692067249215]
eigenvector :: Matrix Double -> Vector Double
eigenvector = fromList . (map realPart) . head . toLists . trans . snd . eig

-- | Calculates the reduced matrix of a given matrix (by reducing the given matrix), e.g.
-- 
-- > reduceMatrix (fromLists [[0.77143,-0.25714],[-0.42245,0.14082]])
--
-- returns
--
-- > (2><2)[ -0.14081563757848092,-0.25714,-0.42245,-0.7714256375784809]
reduceMatrix :: Matrix Double -> Matrix Double
reduceMatrix a = subtractMatrix a (scalarMultiplication (eigenvalue a) (identityMatrix $rows a))

-- | Zipps a matrix col by col
--
-- > zipAllWith sum [[1,2,3],[1,2,3],[1,2,3]] == [3,6,9]
zipAllWith :: (RawVector a -> b) -> RawMatrix a -> RawVector b
zipAllWith _ []  = []
zipAllWith f xss = map f . transpose $ xss