How to extend a matrix in Haskell

Question

I've been trying to create a function that returns a given matrix and adds one row & column whose digits are the sum of the numbers in the row/column next to the digits the sum of the numbers in the column...

E.g.

(1 2
 3 4)

to

(1 2 3
 3 4 7
4 6 10)

My first thought was something like

> extend1 :: Int -> Array (Int,Int) Int
> extend1 ((a,b),(c,d))  = n where
>                 n = array ((a,b),(c,d),(n,n))

because I have to use "array". Now I don't know how to go on. Is it okay or totally wrong?

Answer 1

Using hmatrix:

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util(col,row)

m = (2><2) [1 , 2
           ,3 , 4] :: Matrix Double


extend m = fromBlocks [[m ,sr]
                      ,[sc, s]]
  where
    sr = col $ map sumElements (toRows m)
    sc = row $ map sumElements (toColumns m)
    s  = scalar (sumElements sr)


main = do
    print m
    print $ extend m


(2><2)
 [ 1.0, 2.0
 , 3.0, 4.0 ]
(3><3)
 [ 1.0, 2.0,  3.0
 , 3.0, 4.0,  7.0
 , 4.0, 6.0, 10.0 ]

Answer 2

First of all, you should understand how to interface with haskell arrays. The Array datatype is found in Data.Array, so for full details, look at the docs for that module.

Note that I omit the Ix i constraint you would find on all of these functions. It isn't really important for this situation.

bounds :: Array i e -> (i, i): This function returns the min and max indices of the array. For a 1D array, these are just numbers. For a 2D array, it the the top left and bottom right corners (for a matrix).

array :: (i, i) -> [(i, e)] -> Array i e: This function creates an array from a min/max pair for the bounds, and a list of associations; that is, a map from indices to values. Your intial example can be written as array ((0,0),(1,1)) [((0,0),1),((0,1),2),((1,0),3),((1,1),4)]

assocs :: Array i e -> [(i, e)]: This is the 'opposite' of array. So, arr == array (bounds arr) (assocs arr).

Now the function:

extendArray :: Num e => Array (Int, Int) e -> Array (Int, Int) e
extendArray arr = 
    let
      arr' = assocs arr
      ((xMin, yMin), (xMax, yMax)) = bounds arr
      newCol = [ ((n, yMax + 1) , sum [ v | ((x,_),v) <- arr', x == n] ) | n <- [xMin .. xMax]]
      newRow = [ ((xMax + 1, n) , sum [ v | ((_,y),v) <- arr', y == n] ) | n <- [yMin .. yMax]]
      newCorner = [((xMax + 1, yMax + 1), sum $ map snd arr')]
      newArr = array ((xMin, yMin), (xMax + 1, yMax + 1)) (arr' ++ newCol ++ newRow ++ newCorner)

    in newArr

Lets break down what is happening here. The first two lines in the let statement should be self-explanatory by now. The third line, newCol, is where the magic happens. It uses list comprehensions so if you are unfamiliar with them, see here.

n <- [xMin .. xMax]: This part gets all the possible x-values.

[ v | ((x,_),v) <- arr', x == n]: This reads as, for all values in the list arr' such that the x-coordinate of the index is equal to n, return the value at that coordinate.

((n, yMax + 1) , sum ... Create a new index, whose column is 1 plus the maximum of the old array (so 1 column to the right of the old matrix), whose row is n and whose value is the sum of the previous line (ie - the sum of all the values in row n).

newRow works the same way, but rows and columns reversed.

newCorner gets all the values in the matrix, computes their sum, and places this value at the appropriate index in the new array.

newArr simply combines all the steps into a new array. Obviously the bounds will grow by one in each direction. The associations of the new array are going to include all the old associations, as well as the new ones we have computed.

So:

>let x = array ((0,0),(1,1)) [((0,0),1),((0,1),2),((1,0),3),((1,1),4)]
>x
array ((0,0),(1,1)) [((0,0),1),((0,1),2),((1,0),3),((1,1),4)]
>extendArray x
array ((0,0),(2,2)) [((0,0),1),((0,1),2),((0,2),3),((1,0),3),((1,1),4),((1,2),7),((2,0),4),((2,1),6),((2,2),10)]

Note that this implementation isn't efficient (it is however simple). In newRow, we traverse the entire matrix xMax - xMin times (once for each n value). Since assocs always returns elements in the same order (rows from left to right, columns from up to down) it would be better to split the list arr' into lists each with length yMax - yMin; this would get us a list of the rows. But I will leave this optimization to you.

Answer 3

Using your representation of tuples and without much an eye for generalities, here's a solution to your problem.

type Square2 a = ((a,a), (a,a))
type Square3 a = ((a,a,a), (a,a,a), (a,a,a))

extend1 :: Num a => Square2 a -> Square3 a
extend1 ((a,b), (c,d)) = 
  ( (a,  b,  ab  )
  , (c,  d,  cd  )
  , (ac, bd, abcd) ) 

  where
    ab   = a  + b
    cd   = c  + d
    ac   = a  + c
    bd   = b  + d
    abcd = ab + cd

Note how I use variable definitions culled both from the pattern of the input and the definitions in the where clause. Also notice how I destruct and consume the input while building an entirely new output---some of the elements are reused, but the structure itself is destroyed. Finally note how I chose the type inside of the tuples to be constant and constrained by the Num typeclass which allows me use of (+) for addition.

A similar function, generalized and using Arrays has a type like

extend1A :: Num a => Array (Int,Int) a -> Array (Int, Int) a

where we've lost the ability to know the exact size of the Arrays but it's indicated that our function takes one Array of some size and returns another Array that is indexed identically and containing the same Num-constrained type.

The array function takes as its first argument a set of bounds. We need to update those based on the bounds of the input array

extend1A ary0 = array bounds1 ... where
  ((minX, minY), (maxX, maxY)) = bounds ary0
  (maxX1, maxY1) = (succ maxX, succ maxY)
  bounds1 = ((minX, minY, (maxX1, maxY1))

Then array takes a list of "assocs" as its second argument. We can treat any Array ix a as being equivalent (roughly) to a list of assocs [(ix, a)] which list values stating "at the index ix the value is a". To do this we must know the bounds of the ix type which we managed just previously.

To update an array using the information from an old one, we can modify the assocs of the old array to include new information. Concretely, this means that extend1A looks a bit like

extend1A ary0 = array bounds1 (priorInformation ++ extraInformation) where
  priorInformation = assocs ary0
  extraInformation = ...
  ((minX, minY), (maxX, maxY)) = bounds ary0
  (maxX1, maxY1) = (succ maxX, succ maxY)
  bounds1 = ((minX, minY, (maxX1, maxY1))

If extraInformation were empty ([]) then extend1A ary would be equal to ary on all indicies in range of the input ary and undefined on all outside of its range. We need to fill in extraInformation with the summation information.

extend1A ary0 = array bounds1 (priorInformation ++ extraInformation) where
  priorInformation = assocs ary0 
  extraInformation = xExtension ++ yExtension ++ totalExtension
  xExtension = ...
  yExtension = ...
  totalExtension = ...
  ((minX, minY), (maxX, maxY)) = bounds ary0
  (maxX1, maxY1) = (succ maxX, succ maxY) 
  bounds1 = ((minX, minY, (maxX1, maxY1))

If we think of extending the array in three pieces, the xExtension marked by ab and cd in extend1, the yExtension marked by ac and bd in extend1 and the totalExtension marked by abcd in extend1 we can then compute each part in turn.

totalExtension is easiest--it's just the sum of the "value component" of each (i,a) pair in priorInformation. It could also be the sum of the "value component" of either xExtension or yExtension, but in order to be as obviously correct as possible, we'll pick the first choice and install that at the lower-right corner.

extend1A ary0 = array bounds1 (priorInformation ++ extraInformation) where
  priorInformation = assocs ary0
  extraInformation = xExtension ++ yExtension ++ totalExtension
  sumValues asscs = sum (map snd asscs)
  xExtension = ...
  yExtension = ...
  totalExtension = [((maxX1, maxY1), sumValues priorInformation)]
  ((minX, minY), (maxX, maxY)) = bounds ary0
  (maxX1, maxY1) = (succ maxX, succ maxY)
  bounds1 = ((minX, minY), (maxX1, maxY1))

Note that we can use the where clause to define new functions like sumValues which will show up over and over.

Then we can compute the extensions as list comprehensions over the priorInformation. We need to collect a particular kind of sum over the old assocs---summing over all the values where one index is fixed.

xExtension = [( (maxX1, yix)
              , sumValues (filter (\((_, j), _) -> j == yix) priorInformation)
              )
             | yix <- [minY .. maxY]
             ]

yExtension = [( (xix, maxY1)
              , sumValues (filter (\((i, _), _) -> i == xix) priorInformation)
              )
             | xix <- [minX .. maxX]
             ]

And then we're done. It's not efficient, but all the pieces work together.

extend1A ary0 = array bounds1 (priorInformation ++ extraInformation) where
  priorInformation = assocs ary0
  extraInformation = xExtension ++ yExtension ++ totalExtension
  xExtension = [( (maxX1, yix)
                , sumValues (filter (\((_, j), _) -> j == yix) priorInformation)
                )
               | yix <- [minY .. maxY]
               ]
  yExtension = [( (xix, maxY1)
                , sumValues (filter (\((i, _), _) -> i == xix) priorInformation)
                )
               | xix <- [minX .. maxX]
               ]
  totalExtension = [((maxX1, maxY1), sum xExtension)]
  ((minX, minY), (maxX, maxY)) = bounds ary0
  (maxX1, maxY1) = (succ maxX, succ maxY)
  bounds1 = ((minX, minY), (maxX1, maxY1))