--
-- Copyright (c) 2009-2010, ERICSSON AB All rights reserved.
-- 
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are met:
-- 
--     * Redistributions of source code must retain the above copyright notice,
--       this list of conditions and the following disclaimer.
--     * Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in the
--       documentation and/or other materials provided with the distribution.
--     * Neither the name of the ERICSSON AB nor the names of its contributors
--       may be used to endorse or promote products derived from this software
--       without specific prior written permission.
-- 
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
-- BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
-- OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
-- THE POSSIBILITY OF SUCH DAMAGE.
--

-- | Operations on matrices (doubly-nested parallel vectors). All operations in
-- this module assume rectangular matrices.

module Feldspar.Matrix where



import qualified Prelude as P

import Feldspar.Prelude
import Feldspar.Utils
import Feldspar.Core
import Feldspar.Vector



type Matrix a = Vector (Vector (Data a))



-- | Converts a matrix to a core array.
freezeMatrix :: Storable a => Matrix a -> Data [[a]]
freezeMatrix = freezeVector . map freezeVector

-- | Converts a core array to a matrix. The first length argument is the number
-- of rows (outer vector), and the second argument is the number of columns
-- (inner argument).
unfreezeMatrix :: Storable a => Data Length -> Data Length -> Data [[a]] -> Matrix a
unfreezeMatrix y x = map (unfreezeVector x) . (unfreezeVector y)

-- | Constructs a matrix. The elements are stored in a core array.
matrix :: Storable a => [[a]] -> Matrix a
matrix as
    | allEqual xs = unfreezeMatrix y x (value as)
    | otherwise   = error "matrix: Not rectangular"
  where
    xs = P.map P.length as
    y  = value $ P.length as
    x  = value $ P.head (xs P.++ [0])



-- | Constructing a matrix from an index function.
--
-- @indexedMat m n ixf@:
--
--   * @m@ is the number of rows.
--
--   * @n@ is the number of columns.
--
--   * @ifx@ is a function mapping indexes to elements (first argument is row
--     index; second argument is column index).
indexedMat ::
    Data Int -> Data Int -> (Data Int -> Data Int -> Data a) -> Matrix a

indexedMat m n idx = indexed m $ \k -> indexed n $ \l -> idx k l



-- | Transpose of a matrix
transpose :: Matrix a -> Matrix a
transpose a = indexedMat (length $ head a) (length a) $ \y x -> a ! x ! y
  -- XXX This assumes that (head a) can be used even if a is empty. Might this
  --     violate size constraints on the index?
  --     See the conditional in 'flatten'.



-- | Concatenates the rows of a matrix.
flatten :: Matrix a -> Vector (Data a)
flatten matr = Indexed (m*n) ixf
  where
    m = length matr
    n = (m==0) ? (0, length (head matr))

    ixf i = matr ! y ! x
      where
        y = i `div` n
        x = i `mod` n
  -- XXX Should use "linear indexing"



-- | The diagonal vector of a square matrix. It happens to work if the number of
-- rows is less than the number of columns, but not the other way around (this
-- would require some overhead).
diagonal :: Matrix a -> Vector (Data a)
diagonal m = zipWith (!) m (0 ... (length m - 1))



distributeL :: (a -> b -> c) -> a -> Vector b -> Vector c
distributeL f = map . f

distributeR :: (a -> b -> c) -> Vector a -> b -> Vector c
distributeR = flip . distributeL . flip



class Mul a b
  where
    type Prod a b

    -- | General multiplication operator
    (**) :: a -> b -> Prod a b
      -- XXX This symbol should probably be used for exponentiation instead.

instance Numeric a => Mul (Data a) (Data a)
  where
    type Prod (Data a) (Data a) = Data a
    (**) = (*)

instance Numeric a => Mul (Data a) (DVector a)
  where
    type Prod (Data a) (DVector a) = DVector a
    (**) = distributeL (**)

instance Numeric a => Mul (DVector a) (Data a)
  where
    type Prod (DVector a) (Data a) = DVector a
    (**) = distributeR (**)

instance Numeric a => Mul (Data a) (Matrix a)
  where
    type Prod (Data a) (Matrix a) = Matrix a
    (**) = distributeL (**)

instance Numeric a => Mul (Matrix a) (Data a)
  where
    type Prod (Matrix a) (Data a) = Matrix a
    (**) = distributeR (**)

instance Numeric a => Mul (DVector a) (DVector a)
  where
    type Prod (DVector a) (DVector a) = Data a
    (**) = scalarProd

instance Numeric a => Mul (DVector a) (Matrix a)
  where
    type Prod (DVector a) (Matrix a) = (DVector a)
    vec ** mat = distributeL (**) vec (transpose mat)

instance Numeric a => Mul (Matrix a) (DVector a)
  where
    type Prod (Matrix a) (DVector a) = (DVector a)
    (**) = distributeR (**)

instance Numeric a => Mul (Matrix a) (Matrix a)
  where
    type Prod (Matrix a) (Matrix a) = (Matrix a)
    (**) = distributeR (**)



-- | Matrix multiplication
mulMat :: Numeric a => Matrix a -> Matrix a -> Matrix a
mulMat = (**)

{-# DEPRECATED mul "Please use `mulMat` or `(**)` instead." #-}
-- | Matrix multiplication
mul :: Numeric a => Matrix a -> Matrix a -> Matrix a
mul = (**)



class ElemWise a
  where
    type Elem a

    -- | Operator for general element-wise multiplication
    elemWise :: (Elem a -> Elem a -> Elem a) -> a -> a -> a

instance ElemWise (Data a)
  where
    type Elem (Data a) = Data a
    elemWise = id

instance ElemWise (DVector a)
  where
    type Elem (DVector a) = Data a
    elemWise = zipWith

instance ElemWise (Matrix a)
  where
    type Elem (Matrix a) = Data a
    elemWise = elemWise . elemWise

(.+) :: (ElemWise a, Numeric (Elem a)) => a -> a -> a
(.+) = elemWise (+)

(.-) :: (ElemWise a, Numeric (Elem a)) => a -> a -> a
(.-) = elemWise (-)

(.*) :: (ElemWise a, Numeric (Elem a)) => a -> a -> a
(.*) = elemWise (*)