diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,3 @@ +# Changelog for matrix-lens + +## Unreleased changes diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Author name here (c) 2020 + +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 Author name here nor the names of other + 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 +OWNER 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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,330 @@ +# matrix-lens + +Optics for the `matrix` package. + +![matrix lenses](https://raw.githubusercontent.com/interosinc/matrix-lens/master/matrix-lens.jpg) + +> That’s how it is with people. Nobody cares how it works as long as it works. +> +> - Councillor Hamann + +Sparked by this [reddit post](https://old.reddit.com/r/haskell/comments/gazovx/monthly_hask_anything_may_2020/fqtk9oh/). + +## Examples + +The examples below make use of the following three matrices: + +``` haskell +exampleInt :: Matrix Int +exampleInt = Matrix.fromLists + [ [1, 2, 3] + , [4, 5, 6] + , [7, 8, 9] + ] + +exampleInvertible :: Matrix (Ratio Int) +exampleInvertible = Matrix.fromLists + [ [ -3, 1 ] + , [ 5, 0 ] + ] + +exampleNotSquare :: Matrix Int +exampleNotSquare = Matrix.fromLists + [ [10, 20, 30] + , [40, 50, 60] + , [70, 80, 90] + , [100, 110, 120] + ] +``` + +Accessing individual elements: + +``` haskell +λ> exampleNotSquare ^. elemAt (4, 3) +120 +``` + +``` haskell +λ> exampleInt & elemAt (2, 2) *~ 10 +┌ ┐ +│ 1 2 3 │ +│ 4 50 6 │ +│ 7 8 9 │ +└ ┘ +``` + +Accessing individual columns: + +``` haskell +λ> exampleInt ^. col 1 +[1,4,7] +``` + +``` haskell +λ> exampleInt & col 2 . each *~ 10 +┌ ┐ +│ 1 20 3 │ +│ 4 50 6 │ +│ 7 80 9 │ +└ ┘ +``` + +Accessing individual rows: + +``` haskell +λ> exampleInt ^. row 1 +[1,2,3] +``` + +``` haskell +λ> exampleInt & row 2 . each *~ 100 +┌ ┐ +│ 1 2 3 │ +│ 400 500 600 │ +│ 7 8 9 │ +└ ┘ +``` + +Manipulating all columns as a list: + +``` haskell +λ> exampleInt ^. cols +[[1,4,7],[2,5,8],[3,6,9]] +``` + +``` haskell +λ> exampleInt & cols %~ reverse +┌ ┐ +│ 3 2 1 │ +│ 6 5 4 │ +│ 9 8 7 │ +└ ┘ +``` + +Accessing all rows as a list: + +``` haskell +λ> exampleInt ^. rows +[[1,2,3],[4,5,6],[7,8,9]] +``` + +``` haskell +λ> exampleInt & rows %~ map reverse +┌ ┐ +│ 3 2 1 │ +│ 6 5 4 │ +│ 9 8 7 │ +└ ┘ +``` + +``` haskell +λ> exampleInt & partsOf (dropping 1 (rows . each)) %~ reverse +┌ ┐ +│ 1 2 3 │ +│ 7 8 9 │ +│ 4 5 6 │ +└ ┘ +``` + +In addition to the above there are also `switching` and `sliding` Isos for both +rows and columns which allow you to swap two arbitrary rows or columns or slide +a row or column through the matrix to a different row or column (moving all +intervening rows or columns over in the direction of the source row or column): + +``` haskell +λ> exampleNotSquare ^. switchingRows 1 4 +┌ ┐ +│ 100 110 120 │ +│ 40 50 60 │ +│ 70 80 90 │ +│ 10 20 30 │ +└ ┘ +``` + +``` haskell +λ> exampleNotSquare ^. slidingRows 1 4 +┌ ┐ +│ 40 50 60 │ +│ 70 80 90 │ +│ 100 110 120 │ +│ 10 20 30 │ +└ ┘ +``` + +..and similary for `switchingCols` and `switchingRows`. + +An Iso exists for accessing the matrix with a given row scaled: + +``` haskell +λ> exampleInt ^. scalingRow 1 10 +┌ ┐ +│ 10 20 30 │ +│ 4 5 6 │ +│ 7 8 9 │ +└ ┘ + +λ> exampleInt & scalingRow 1 10 . flattened *~ 2 +┌ ┐ +│ -200 -400 -600 │ +│ 8 10 12 │ +│ 14 16 18 │ +└ ┘ +``` + +Any valid sub matrix can be accessed via the `sub` lens: + +``` haskell +λ> exampleNotSquare ^. sub (2, 1) (3, 2) +┌ ┐ +│ 40 50 │ +│ 70 80 │ +└ ┘ + +λ> exampleNotSquare & sub (2, 1) (3, 2) . rows %~ reverse +┌ ┐ +│ 10 20 30 │ +│ 70 80 60 │ +│ 40 50 90 │ +│ 100 110 120 │ +└ ┘ +``` + +The transposition of the matrix can be accessed via the `transposed` Iso: + +``` haskell +λ> exampleInt ^. transposed +┌ ┐ +│ 1 4 7 │ +│ 2 5 8 │ +│ 3 6 9 │ +└ ┘ + +λ> exampleInt & transposed . taking 4 flattened *~ 10 +┌ ┐ +│ 10 20 3 │ +│ 40 5 6 │ +│ 70 8 9 │ +└ ┘ +``` + + + +You can also traverse the `flattened` matrix: + +``` haskell +λ> exampleInt ^.. flattened +[1,2,3,4,5,6,7,8,9] +``` + +which is more useful for making modifications: + +``` haskell +λ> exampleInt & flattened . filtered even *~ 10 +┌ ┐ +│ 1 20 3 │ +│ 40 5 60 │ +│ 7 80 9 │ +└ ┘ +``` + +``` haskell +λ> exampleInt & dropping 4 flattened *~ 10 +┌ ┐ +│ 1 2 3 │ +│ 4 50 60 │ +│ 70 80 90 │ +└ ┘ +``` + +Accessing the diagonal: + +``` haskell +λ> exampleInt ^. diag +[1,5,9] + +λ> exampleInt & diag %~ reverse +┌ ┐ +│ 9 2 3 │ +│ 4 5 6 │ +│ 7 8 1 │ +└ ┘ + +λ> exampleInt & diag . each *~ 10 +┌ ┐ +│ 10 2 3 │ +│ 4 50 6 │ +│ 7 8 90 │ +└ ┘ +``` + +Accessing inverse matrix is possible via the `inverted` optic. Since not all +matrices have inverses `inverted` is a prism: + +``` haskell +λ> exampleInvertible ^? inverted +Just ┌ ┐ +│ 0 % 1 1 % 5 │ +│ 1 % 1 3 % 5 │ +└ ┘ + +λ> exampleInvertible & inverted . flattened *~ 2 +┌ ┐ +│ (-3) % 2 1 % 2 │ +│ 5 % 2 0 % 1 │ +└ ┘ +``` + +Minor matrices can be accessed by specifying the (r, c) to be removed: + +``` haskell +λ> exampleInt ^. minor (1, 2) +┌ ┐ +│ 4 6 │ +│ 7 9 │ +└ ┘ + +λ> exampleInt & minor (1, 2) . flattened *~ 10 +┌ ┐ +│ 1 2 3 │ +│ 40 5 60 │ +│ 70 8 90 │ +└ ┘ +``` + +An Iso exists for accessing a scaled version of a matrix: + +``` haskell +λ> exampleInt ^. scaled 10 +┌ ┐ +│ 10 20 30 │ +│ 40 50 60 │ +│ 70 80 90 │ +└ ┘ + +λ> exampleInt & minor (1, 1) . scaled 10 . flattened +~ 1 +┌ ┐ +│ 1 2 3 │ +│ 4 -510 -610 │ +│ 7 -810 -910 │ +└ ┘ +``` + +Getters for the matrix determinant and size are also provided: + +``` haskell +λ> exampleInt ^. determinant +Just 0 + +λ> exampleInvertible ^. determinant +Just ((-5) % 1) + +λ> exampleNotSquare ^. determinant +Nothing + +λ> exampleInt ^. size +(3,3) +λ> exampleInvertible ^. size +(2,2) +λ> exampleNotSquare ^. size +(4,3) +``` diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/matrix-lens.cabal b/matrix-lens.cabal new file mode 100644 --- /dev/null +++ b/matrix-lens.cabal @@ -0,0 +1,79 @@ +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.33.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: 853a69950c400cb9cafc73a36a7bfd908a1c2890c5a47196defaeb91ba08e83f + +name: matrix-lens +version: 0.1.0.0 +synopsis: Optics for the "matrix" package +description: See README at +category: Math +homepage: https://github.com/interosinc/matrix-lens#readme +bug-reports: https://github.com/interosinc/matrix-lens/issues +author: Interos, Inc. +maintainer: jevans@interos.ai +copyright: 2020 Interos, Inc. +license: BSD3 +license-file: LICENSE +build-type: Simple +extra-source-files: + README.md + ChangeLog.md + +source-repository head + type: git + location: https://github.com/interosinc/matrix-lens + +flag developer + description: Developer mode -- stricter handling of compiler warnings. + manual: True + default: False + +library + exposed-modules: + Data.Matrix.Lens + Data.Matrix.Lens.Internal + Data.Matrix.Lens.Internal.Warnings + other-modules: + Paths_matrix_lens + hs-source-dirs: + src + build-depends: + base >=4.7 && <5 + , lens >=4.19.2 && <4.20 + , matrix >=0.3.6 && <0.4 + , vector >=0.12.1 && <0.13 + if flag(developer) + ghc-options: -Weverything -Werror -Wno-all-missed-specialisations -Wno-missed-specialisations -Wno-missing-deriving-strategies -Wno-missing-exported-signatures -Wno-missing-import-lists -Wno-missing-local-signatures -Wno-monomorphism-restriction -Wno-safe -Wno-unsafe + else + ghc-options: -Weverything -Wno-all-missed-specialisations -Wno-missed-specialisations -Wno-missing-deriving-strategies -Wno-missing-exported-signatures -Wno-missing-import-lists -Wno-missing-local-signatures -Wno-monomorphism-restriction -Wno-safe -Wno-unsafe + default-language: Haskell2010 + +test-suite matrix-lens-test + type: exitcode-stdio-1.0 + main-is: Driver.hs + other-modules: + MatrixLensTest + Paths_matrix_lens + hs-source-dirs: + test + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.7 && <5 + , hedgehog + , lens >=4.19.2 && <4.20 + , matrix >=0.3.6 && <0.4 + , matrix-lens + , tasty + , tasty-discover + , tasty-hedgehog + , tasty-hspec + , vector >=0.12.1 && <0.13 + if flag(developer) + ghc-options: -Weverything -Werror -Wno-all-missed-specialisations -Wno-missed-specialisations -Wno-missing-deriving-strategies -Wno-missing-exported-signatures -Wno-missing-import-lists -Wno-missing-local-signatures -Wno-monomorphism-restriction -Wno-safe -Wno-unsafe + else + ghc-options: -Weverything -Wno-all-missed-specialisations -Wno-missed-specialisations -Wno-missing-deriving-strategies -Wno-missing-exported-signatures -Wno-missing-import-lists -Wno-missing-local-signatures -Wno-monomorphism-restriction -Wno-safe -Wno-unsafe + default-language: Haskell2010 diff --git a/src/Data/Matrix/Lens.hs b/src/Data/Matrix/Lens.hs new file mode 100644 --- /dev/null +++ b/src/Data/Matrix/Lens.hs @@ -0,0 +1,85 @@ +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Data.Matrix.Lens + ( flattened + , col + , cols + , determinant + , diag + , elemAt + , inverted + , isSquare + , minor + , row + , rows + , scaled + , scalingRow + , sub + , size + , slidingCols + , slidingRows + , switchingCols + , switchingRows + , transposed + ) where + +import Prelude + +import Control.Lens hiding ( set ) +import Data.Bifunctor ( first ) +import qualified Data.Foldable as F ( toList ) +import qualified Data.List as L +import Data.Matrix +import Data.Matrix.Lens.Internal as X ( col + , elemAt + , isSquare + , minor + , row + , rows + , slidingCols + , slidingRows + , switchingCols + , switchingRows + ) +import Data.Matrix.Lens.Internal.Warnings ( determinant + , size + ) +import Data.Maybe ( fromMaybe ) +import qualified Data.Vector as V + +transposed :: Iso' (Matrix a) (Matrix a) +transposed = iso transpose transpose + +scaled :: Num a => a -> Iso' (Matrix a) (Matrix a) +scaled n = iso (scaleMatrix n) (scaleMatrix . negate $ n) + +cols :: Lens' (Matrix a) [[a]] +cols = lens (L.transpose . toLists) (const (fromLists . L.transpose)) + +scalingRow :: Num a => Int -> a -> Iso' (Matrix a) (Matrix a) +scalingRow r n = iso (scaleRow n r) (scaleRow (negate n) r) + +flattened :: Traversal' (Matrix a) a +flattened = rows . each . each + +sub :: (Int, Int) -> (Int, Int) -> Lens' (Matrix a) (Matrix a) +sub (r1, c1) (r2, c2) = lens (submatrix r1 r2 c1 c2) (setSubMatrix (r1, c1)) + where + setSubMatrix (r, c) dst src = foldr f dst indexedRows + where + indexedRows = zip [r..] . map (zip [c..]) . toLists $ src + f (r', indexedCols) dst' = foldr (g r') dst' indexedCols + g r' (c', x) dst' = fromMaybe dst' $ safeSet x (r', c') dst' + +inverted :: (Eq a, Fractional a) => Prism' (Matrix a) (Matrix a) +inverted = flip prism (\x -> first (const x) . inverse $ x) $ \x -> case inverse x of + Left _ -> x + Right y -> y + +diag :: Lens' (Matrix a) [a] +diag = lens (V.toList . getDiag) (\m -> setDiag m . V.fromList) + where + setDiag m = foldr f m . zip [1..] . F.toList + where + f (n, x) m' = m' & elemAt (n, n) .~ x diff --git a/src/Data/Matrix/Lens/Internal.hs b/src/Data/Matrix/Lens/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Matrix/Lens/Internal.hs @@ -0,0 +1,120 @@ +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# OPTIONS_GHC -Wno-missing-export-lists #-} + +module Data.Matrix.Lens.Internal where + +import Prelude + +import Control.Lens +import Data.Matrix +import Data.Maybe ( fromMaybe ) +import Data.Vector ( Vector ) +import qualified Data.Vector as V + +isSquare :: Matrix a -> Bool +isSquare = uncurry (==) . getSize + +getSize :: Matrix a -> (Int, Int) +getSize m = (nrows m, ncols m) + +elemAt :: (Int, Int) -> Lens' (Matrix a) a +elemAt (i, j) = lens (getElem i j) (\m x -> setElem x (i, j) m) + +row :: Int -> Lens' (Matrix a) [a] +row r = lens (V.toList . getRow r) (\m -> setRow r m . V.fromList) + +col :: Int -> Lens' (Matrix a) [a] +col c = lens (V.toList . getCol c) (\m -> setCol c m . V.fromList) + +rows :: Lens' (Matrix a) [[a]] +rows = lens toLists (const fromLists) + +switchingRows :: Int -> Int -> Iso' (Matrix a) (Matrix a) +switchingRows r1 r2 = iso (switchRows r1 r2) (switchRows r2 r1) + +switchingCols :: Int -> Int -> Iso' (Matrix a) (Matrix a) +switchingCols c1 c2 = iso (switchCols c1 c2) (switchCols c2 c1) + +slidingRows :: Int -> Int -> Iso' (Matrix a) (Matrix a) +slidingRows r1 r2 = iso (slideRows r1 r2) (slideRows r2 r1) + +slidingCols :: Int -> Int -> Iso' (Matrix a) (Matrix a) +slidingCols c1 c2 = iso (slideCols c1 c2) (slideCols c2 c1) + +minor :: (Int, Int) -> Lens' (Matrix a) (Matrix a) +minor (r, c) = lens (minorMatrix r c) (setMinorMatrix (r, c)) + +-- ================================================================ -- + +getDeterminant :: Num a => Matrix a -> Maybe a +getDeterminant m + | getSize m == (2, 2) = Just . twoByTwoDet $ m + | otherwise = laplace m + +twoByTwoDet :: Num a => Matrix a -> a +twoByTwoDet m = case map ((m ^.) . elemAt) coords of + [a, b, c, d] -> a * d - b * c + _other -> error "unpossible! (2)" + where + coords = + [ (1, 1) + , (1, 2) + , (2, 1) + , (2, 2) + ] + +laplace :: Num a => Matrix a -> Maybe a +laplace m + | not . isSquare $ m = Nothing + | otherwise = Just . sum . zipWith (*) (cycle [1, -1]) . map f $ [1..ncols m] + where + r = 1 + f c = ((e *) <$> getDeterminant (m ^. minor pair)) + & fromMaybe (error "unpossible! (3)") + where + e = m ^. elemAt pair + pair = (r, c) + +-- ================================================================ -- + +setRow :: Int -> Matrix a -> Vector a -> Matrix a +setRow r m v = foldr (\(c, x) -> setElem x (r, c)) m $ + zip [1..] (V.toList v) + +setCol :: Int -> Matrix a -> Vector a -> Matrix a +setCol c m v = foldr (\(r, x) -> setElem x (r, c)) m $ + zip [1..] (V.toList v) + +setMinorMatrix :: forall a. (Int, Int) -> Matrix a -> Matrix a -> Matrix a +setMinorMatrix (r, c) dst src = sequenceA inserted + & fromMaybe (error "unpossible! (1)") + where + inserted = foldr copyCol m' indexedCol + where + m' = foldr copyRow adjusted indexedRow + + indexedCol = zip [1..] $ dst ^. col c + indexedRow = zip [1..] $ dst ^. row r + + copyRow (c', x) = elemAt (r, c') ?~ x + copyCol (r', x) = elemAt (r', c) ?~ x + + adjusted = injected ^. slidingRows (nrows mm + 1) r + . slidingCols (ncols mm + 1) c + + injected = extendTo Nothing (nrows mm + 1) (ncols mm + 1) mm + + mm = Just <$> src + +slideRows :: Int -> Int -> Matrix a -> Matrix a +slideRows s d m + | s > d = slideRows (s - 1) d $ m ^. switchingRows s (s - 1) + | s < d = slideRows (s + 1) d $ m ^. switchingRows s (s + 1) + | otherwise = m + +slideCols :: Int -> Int -> Matrix a -> Matrix a +slideCols s d m + | s > d = slideCols (s - 1) d $ m ^. switchingCols s (s - 1) + | s < d = slideCols (s + 1) d $ m ^. switchingCols s (s + 1) + | otherwise = m diff --git a/src/Data/Matrix/Lens/Internal/Warnings.hs b/src/Data/Matrix/Lens/Internal/Warnings.hs new file mode 100644 --- /dev/null +++ b/src/Data/Matrix/Lens/Internal/Warnings.hs @@ -0,0 +1,28 @@ +-- These functions are in a separate file so we can cover only these two +-- functions with the -Wno-redundant-constraints pragma and leave +-- -Wredundant-constraints enabled for everything else. If anyone knows a +-- better way to handle this, please let me know. + +{-# OPTIONS_GHC -Wno-redundant-constraints #-} + +module Data.Matrix.Lens.Internal.Warnings + ( determinant + , getSize + , size + ) where + +import Prelude + +import Control.Lens ( Getter + , to + ) +import Data.Matrix ( Matrix ) +import Data.Matrix.Lens.Internal ( getDeterminant + , getSize + ) + +determinant :: Num a => Getter (Matrix a) (Maybe a) +determinant = to getDeterminant + +size :: Getter (Matrix a) (Int, Int) +size = to getSize diff --git a/test/Driver.hs b/test/Driver.hs new file mode 100644 --- /dev/null +++ b/test/Driver.hs @@ -0,0 +1,1 @@ +{-# OPTIONS_GHC -F -pgmF tasty-discover #-} diff --git a/test/MatrixLensTest.hs b/test/MatrixLensTest.hs new file mode 100644 --- /dev/null +++ b/test/MatrixLensTest.hs @@ -0,0 +1,431 @@ +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module MatrixLensTest + ( hprop_diagIsLesserOfRC + , hprop_invertedHasSameDimensions + , hprop_nonSquareMatricesHaveNoDeterminants + , hprop_squareMatricesHaveDeterminants + , spec_determinant + , spec_diag + , spec_elemAt + , spec_examples + , spec_inverted + , spec_minor + , spec_row + , spec_sub + ) where + +import Prelude + +import Control.Lens ( (%~) + , (&) + , (*~) + , (.~) + , (^.) + , (^?) + , Lens' + , each + , partsOf + , set + , view + ) +import Control.Monad ( guard + , replicateM + ) +import Data.Foldable ( traverse_ ) +import Data.Matrix ( Matrix ) +import qualified Data.Matrix as Matrix +import Data.Matrix.Lens +import Data.Maybe ( isJust + , isNothing + ) +import Data.Ratio ( (%) + , Ratio + ) +import Hedgehog ( (===) + , Gen + , MonadGen + , Property + , assert + , forAll + , property + , withDiscards + ) +import qualified Hedgehog.Gen as Gen +import qualified Hedgehog.Range as Range +import Test.Tasty.Hspec + +spec_elemAt :: Spec +spec_elemAt = do + let m = exampleInt + + context "views the appropriate locations" $ do + + let testView (pair, value) = let (label, p) = setup pair in + it ("at " <> label) $ + m ^. elemAt p `shouldBe` value + + traverse_ testView + [ ((1, 1), 1) + , ((1, 2), 2) + , ((2, 1), 4) + , ((2, 2), 5) + ] + + context "sets the appropriate locations" $ do + let testSet (pair, expected) = let (label, p) = setup pair in + it ("at " <> label) $ + (m & elemAt p .~ 99) `shouldBe` Matrix.fromLists expected + + traverse_ testSet + [ ((1, 1), [ [99, 2, 3] + , [4, 5, 6] + , [7, 8, 9] + ]) + , ((1, 2), [ [1, 99, 3] + , [4, 5, 6] + , [7, 8, 9] + ]) + , ((2, 1), [ [1, 2, 3] + , [99, 5, 6] + , [7, 8, 9] + ]) + , ((2, 2), [ [1, 2, 3] + , [4, 99, 6] + , [7, 8, 9] + ]) + ] + +spec_row :: Spec +spec_row = do + let m = exampleInt + + it "views the appropriate rows" $ + m ^. row 1 `shouldBe` [1, 2, 3] + +spec_sub :: Spec +spec_sub = do + let m = exampleInt + + it "sets the appropriate locations" $ + (m & sub (2, 2) (3, 3) .~ m) `shouldBeMatrix` + [ [1, 2, 3] + , [4, 1, 2] + , [7, 4, 5] + ] + + it "modifies the appropriate locations" $ + (m & sub (2, 2) (3, 3) %~ Matrix.transpose) `shouldBeMatrix` + [ [1, 2, 3] + , [4, 5, 8] + , [7, 6, 9] + ] + +spec_minor :: Spec +spec_minor = do + + context "reads the appropriate locations" $ do + + let testView (pair, value) = let (label, p) = setup pair in + it ("at " <> label) $ + exampleInt ^. minor p `shouldBeMatrix` value + + traverse_ testView + [ ((1, 1), [ [5, 6] + , [8, 9] + ]) + , ((2, 2), [ [1, 3] + , [7, 9] + ]) + , ((1, 2), [ [4, 6] + , [7, 9] + ]) + ] + + context "sets the appropriate locations" $ do + + let testSet (pair, value) = let (label, p) = setup pair in + it ("at " <> label) $ + (exampleInt & minor p %~ Matrix.transpose) `shouldBeMatrix` value + + traverse_ testSet + [ ((1, 1), [ [ 1, 2, 3 ] + , [ 4, 5, 8 ] + , [ 7, 6, 9 ] + ]) + , ((2, 2), [ [ 1, 2, 7 ] + , [ 4, 5, 6 ] + , [ 3, 8, 9 ] + ]) + , ((1, 2), [ [ 1, 2, 3 ] + , [ 4, 5, 7 ] + , [ 6, 8, 9 ] + ]) + ] + +spec_inverted :: Spec +spec_inverted = do + it "inverts an invertible matrix" $ + exampleInvertible ^? inverted `shouldBeJustMatrix` + [ [ 0 % 1, 1 % 5 ] + , [ 1 % 1, 3 % 5 ] + ] + + it "roundtrips" $ + exampleInvertible ^? inverted . inverted `shouldBe` Just + exampleInvertible + + it "modifies correctly" $ + (exampleInvertible & inverted . elemAt (1, 1) .~ 5 % 2) `shouldBeMatrix` + [ [ 6 % 13, (-2) % 13 ] + , [ (-10) % 13, 25 % 13 ] + ] + +hprop_invertedHasSameDimensions :: Property +hprop_invertedHasSameDimensions = withDiscards 200 . property $ do + m <- forAll genInvertibleMatrix + m ^? inverted . size === Just (m ^. size) + +spec_diag :: Spec +spec_diag = do + context "given a square matrix" $ do + let m = exampleInt + + it "reads the right values" $ + m ^. diag `shouldBe` [1, 5, 9] + + it "writes the right values" $ + (m & diag .~ [20, 60, 100]) `shouldBeMatrix` + [ [ 20, 2, 3 ] + , [ 4, 60, 6 ] + , [ 7, 8, 100 ] + ] + + context "given a non-square matrix" $ do + let m = exampleNotSquare + + it "reads the right values" $ + m ^. diag `shouldBe` [10, 50, 90] + + + it "writes the right values" $ + (m & diag .~ [1, 2, 3]) `shouldBeMatrix` + [ [ 1, 20, 30] + , [ 40, 2, 60] + , [ 70, 80, 3] + , [100, 110, 120] + ] + +hprop_diagIsLesserOfRC :: Property +hprop_diagIsLesserOfRC = property $ do + m <- forAll $ Gen.choice [genSquareMatrix, genNonSquareMatrix] + length (m ^. diag) === min (Matrix.nrows m) (Matrix.ncols m) + +spec_examples :: Spec +spec_examples = do + it "should be able to transpose a minor matrix" $ + (exampleInt & minor (1, 1) %~ Matrix.transpose) `shouldBeMatrix` + [ [1, 2, 3] + , [4, 5, 8] + , [7, 6, 9] + ] + + it "should be able to reverse rows" $ + (exampleInt & rows %~ reverse) `shouldBeMatrix` + [ [7, 8, 9] + , [4, 5, 6] + , [1, 2, 3] + ] + + it "should be able to reverse cols" $ + (exampleInt & cols %~ reverse) `shouldBeMatrix` + [ [3, 2, 1] + , [6, 5, 4] + , [9, 8, 7] + ] + + it "should be able to set a minor matrix to one value" $ + (exampleInt & minor (2, 2) . flattened .~ 1) `shouldBeMatrix` + [ [1, 2, 1] + , [4, 5, 6] + , [1, 8, 1] + ] + + it "should be able to set everything top to bottom" $ + (exampleInt & partsOf flattened .~ [90,80..]) `shouldBeMatrix` + [ [90, 80, 70] + , [60, 50, 40] + , [30, 20, 10] + ] + + it "should be able to reverse all cells" $ + (exampleInt & partsOf flattened %~ reverse) `shouldBeMatrix` + [ [9, 8, 7] + , [6, 5, 4] + , [3, 2, 1] + ] + +spec_determinant :: Spec +spec_determinant = do + let m3x3 = Matrix.fromLists + [ [6, 1, 1 :: Int] + , [4, -2, 5] + , [2, 8, 7] + ] + m2x2 = m3x3 ^. minor (1, 1) + + it "should return Nothing on non-square matrices" $ + exampleNotSquare ^. determinant `shouldBe` Nothing + + it "should work for 2x2 matrices" $ + m2x2 ^. determinant `shouldBe` Just (-54) + + it "should work for > 2x2 square matrices" $ + m3x3 ^. determinant `shouldBe` Just (-306) + +hprop_squareMatricesHaveDeterminants :: Property +hprop_squareMatricesHaveDeterminants = property $ do + m <- forAll genSquareMatrix + assert . isJust $ m ^. determinant + +hprop_nonSquareMatricesHaveNoDeterminants :: Property +hprop_nonSquareMatricesHaveNoDeterminants = property $ do + m <- forAll genNonSquareMatrix + assert . isNothing $ m ^. determinant + +-- ================================================================ -- + +infix 1 `shouldBeMatrix` +shouldBeMatrix :: (Eq a, Show a) => Matrix a -> [[a]] -> Expectation +shouldBeMatrix x y = x `shouldBe` Matrix.fromLists y + +infix 1 `shouldBeJustMatrix` +shouldBeJustMatrix :: (Eq a, Show a) => Maybe (Matrix a) -> [[a]] -> Expectation +shouldBeJustMatrix x y = x `shouldBe` Just (Matrix.fromLists y) + +setup :: (Int, Int) -> (String, (Int, Int)) +setup = (,) =<< show + +-- ================================================================ -- + +genSquareMatrix :: Gen (Matrix Int) +genSquareMatrix = do + sz <- genSize + flip (set (partsOf flattened)) (Matrix.identity sz) <$> genValues (sz, sz) + +genNonSquareMatrix :: Gen (Matrix Int) +genNonSquareMatrix = do + r <- genSize + c <- genSize + guard $ r /= c + let m = Matrix.extendTo 0 r c . Matrix.identity $ min r c + flip (set (partsOf flattened)) m <$> genValues (r, c) + +type MRI = Matrix RI +type RI = Ratio Int + +data ElementaryOp + = InterchangeCols Int Int + | InterchangeRows Int Int + | ScaleRow Int RI + | ScaleCol Int RI + | ScaleAndAddRow Int Int RI + | ScaleAndAddCol Int Int RI + deriving (Eq, Show) + +genInvertibleMatrix :: Gen MRI +genInvertibleMatrix = do + (sz, im) <- (,) <*> Matrix.identity <$> genSize + nOps <- genSize + foldr ($) im <$> replicateM nOps (genOp sz) + where + genOp :: Int -> Gen (MRI -> MRI) + genOp sz = makeFun <$> genEOp sz + + makeFun :: ElementaryOp -> (MRI -> MRI) + makeFun = \case + InterchangeRows r1 r2 -> view (switchingRows r1 r2) + InterchangeCols c1 c2 -> view (switchingCols c1 c2) + ScaleRow r n -> row r . each *~ n + ScaleCol c n -> col c . each *~ n + ScaleAndAddRow r1 r2 n -> scaleAndAdd row r1 r2 n + ScaleAndAddCol c1 c2 n -> scaleAndAdd col c1 c2 n + + scaleAndAdd :: (Int -> Lens' MRI [RI]) -> Int -> Int -> RI -> MRI -> MRI + scaleAndAdd acc a b n m = m & acc a %~ zipWith (+) (map (*n) $ m ^. acc b) + + genEOp :: Int -> Gen ElementaryOp + genEOp n = Gen.choice + [ genIR n + , genIC n + , genSR n + , genSC n + , genAR n + , genAC n + ] + + genIR :: Int -> Gen ElementaryOp + genIR n = do + r1 <- genOneToN n + r2 <- genOneToN n + guard $ r1 /= r2 + pure $ InterchangeRows r1 r2 + + genIC :: Int -> Gen ElementaryOp + genIC n = do + c1 <- genOneToN n + c2 <- genOneToN n + guard $ c1 /= c2 + pure $ InterchangeCols c1 c2 + + genSR :: Int -> Gen ElementaryOp + genSR n = ScaleRow <$> genOneToN n <*> genScale + + genSC :: Int -> Gen ElementaryOp + genSC n = ScaleRow <$> genOneToN n <*> genScale + + genAR :: Int -> Gen ElementaryOp + genAR n = ScaleAndAddRow <$> genOneToN n <*> genOneToN n <*> genScale + + genAC :: Int -> Gen ElementaryOp + genAC n = ScaleAndAddCol <$> genOneToN n <*> genOneToN n <*> genScale + + genOneToN :: Int -> Gen Int + genOneToN n = Gen.int (Range.linearFrom 1 1 n) + + genScale :: Gen RI + genScale = fromIntegral <$> Gen.int (Range.linearFrom 1 1 500) + +genValues :: (MonadGen m, Integral a) => (Int, Int) -> m [a] +genValues (r, c) = Gen.list (Range.singleton $ r * c) genInt + where + genInt = Gen.integral (Range.linearFrom 0 (-100) 100) + +genSize :: MonadGen m => m Int +genSize = Gen.integral (Range.linear 2 10) + +-- ================================================================ -- + +exampleInt :: Matrix Int +exampleInt = Matrix.fromLists + [ [1, 2, 3] + , [4, 5, 6] + , [7, 8, 9] + ] + +exampleInvertible :: Matrix (Ratio Int) +exampleInvertible = Matrix.fromLists + [ [ -3, 1 ] + , [ 5, 0 ] + ] + +exampleNotSquare :: Matrix Int +exampleNotSquare = Matrix.fromLists + [ [10, 20, 30] + , [40, 50, 60] + , [70, 80, 90] + , [100, 110, 120] + ]