Merge pull request #787 from explorable-viz/example-fluid

MatrixUpdate Created

example/Example.purs

@@ -2,10 +2,11 @@ module Example.Example where
 
 import Prelude
 
-import Data.Array ((..))
+-- import Data.Array ((..))
+import Data.List ((:), List(..))
 import Effect (Effect)
 import Effect.Class.Console (log, logShow)
-import Example.Util.DTW (distEuclid, distanceDTWWindow)
+import Example.Util.DTW (costMatrixInit, distEuclid, distanceDTWWindow)
 import Partial.Unsafe (unsafePartial)
 import Util ((×))
 
@@ -16,14 +17,11 @@ main = do
       n = 5
       m = 7
       window = 2
-      nextIndices = do
-         i <- 1 .. n
-         j <- (max 1 (i - window)) .. (min m (i + window))
-         [ (i × j) ]
-   logShow nextIndices
+      initMat = costMatrixInit n m window
+   logShow initMat
    let
-      x = [ 3.0, 1.0, 2.0, 2.0, 1.0 ]
-      y = [ 2.0, 0.0, 0.0, 3.0, 3.0, 1.0, 0.0 ]
+      x = (3.0 : 1.0 : 2.0 : 2.0 : 1.0 : Nil)
+      y = (2.0 : 0.0 : 0.0 : 3.0 : 3.0 : 1.0 : 0.0 : Nil)
       m1 × m2 = unsafePartial $ distanceDTWWindow x y 2 distEuclid
    log "Finished DTW"
    logShow m2

example/Util/DTW.purs

@@ -1,14 +1,15 @@
 module Example.Util.DTW
    ( NumInf(..)
+   , costMatrixInit
    , distEuclid
    , distanceDTWWindow
    ) where
 
 import Prelude
 
-import Data.Array (length, modifyAtIndices, range, (!!), (..))
+import Data.Array (modifyAtIndices, range, (!!), (..))
 import Data.Foldable (foldl)
-import Data.List (List(..))
+import Data.List (List(..), index, length)
 import Data.Maybe (fromJust)
 import Data.Ord (abs)
 import Partial.Unsafe (unsafePartial)
@@ -29,7 +30,7 @@ costMatrixInit rows cols window = mapMatrix init indexMat
    init (_ × 0) = Infty
    init (x × y) = if (abs $ x - y) <= window then FNum 0.0 else Infty
 
-distanceDTWWindow :: Partial => Array Number -> Array Number -> Int -> (Number -> Number -> NumInf) -> Matrix NumInf × List (Int × Int)
+distanceDTWWindow :: Partial => List Number -> List Number -> Int -> (Number -> Number -> NumInf) -> Matrix NumInf × List (Int × Int)
 distanceDTWWindow seq1 seq2 window cost = result × (extractPath priorcells (n × m))
    where
    n = length seq1
@@ -48,7 +49,7 @@ distanceDTWWindow seq1 seq2 window cost = result × (extractPath priorcells (n 
          ijm1 = dists ! i' ! j' - 1
          im1jm1 = dists ! i' - 1 ! j' - 1
          minim × prev = costAndPrevD (i' × j') im1j ijm1 im1jm1
-         costij = cost (seq1 ! i' - 1) (seq2 ! j' - 1) `plus` minim
+         costij = cost (seq1 ? i' - 1) (seq2 ? j' - 1) `plus` minim
       in
          updateAt i' j' dists (const costij) × updateAt i' j' inds (const prev)
 
@@ -92,6 +93,11 @@ unsafeArrayInd arr ind = unsafePartial $ fromJust (arr !! ind)
 
 infixl 5 unsafeArrayInd as !
 
+unsafeListInd :: forall a. List a -> Int -> a
+unsafeListInd list ind = unsafePartial $ fromJust (index list ind)
+
+infixl 5 unsafeListInd as ?
+
 mapMatrix :: forall a b. (a -> b) -> Matrix a -> Matrix b
 mapMatrix f m = map (\row -> map f row) m
 

fluid/example/dtw/cost-matrix.fld

@@ -0,0 +1,11 @@
+let abs x y = if x - y < 0 then negate (x - y) else (x - y);
+    costMatrixInit rows cols window = 
+      let ns = [0 .. rows];
+          ms = [0 .. cols]
+      in 
+        [| let initV = if or (and (n == 1) (m == 1)) (and (abs n m <= window) (not (or (n == 1) (m == 1))))
+                       then FNum 0
+                       else Infty 
+           in initV | (n, m) in (rows, cols) |]
+in costMatrixInit 5 7 1
+

fluid/example/dtw/cost-predicate.fld

@@ -0,0 +1,5 @@
+let abs x y = if x - y < 0 then negate (x - y) else (x - y);
+    initV n m window = if (or (and (n == 0) (m == 0)) (and ((abs n m) <= window) (not (or (n == 0) (m == 0)))))
+                       then FNum 0 
+                       else Infty 
+in initV 1 0 2

fluid/example/dtw/matrix-update.fld

@@ -0,0 +1,2 @@
+let mat = [| (n, m) | (n,m) in (3, 3)|]
+in matrixUpdate mat (1,1) 100

fluid/example/dtw/min-prev.fld

@@ -0,0 +1,10 @@
+let minAndPrev (i, j) im1 jm1 ijm1 = 
+    let minimal = minimum [im1,jm1,ijm1]
+    in 
+    if minimal == im1 
+    then (((i - 1), j), minimal)
+    else if minimal == jm1
+    then ((i, (j - 1)), minimal)
+    else
+        (((i - 1), (j - 1)), ijm1)
+in minAndPrev (4,5) 3 4 5

fluid/example/dtw/next-indices.fld

@@ -0,0 +1,5 @@
+let nextIndices n m window = 
+    [(i, j) | i <- [1 .. n],
+              j <- [(max 1 (i - window)) .. (min m (i + window))]]
+in 
+  nextIndices 5 7 2

src/DataType.purs

@@ -168,4 +168,8 @@ dataTypes = L.fromFoldable
    , dataType "Marker"
         [ "Arrowhead" × 0
         ]
+   , dataType "InfNum"
+        [ "FNum" × 1
+        , "Infty" × 0
+        ]
    ]

src/Primitive/Defs.purs

@@ -66,6 +66,7 @@ primitives = D.fromFoldable
    , "mod" × binaryZero { i: int, o: int, fwd: mod }
    , "quot" × binaryZero { i: int, o: int, fwd: quot }
    , "rem" × binaryZero { i: int, o: int, fwd: rem }
+   , "matrixUpdate" × extern matrixMut
    ]
 
 error_ :: ForeignOp
@@ -96,8 +97,6 @@ debugLog = mkExists $ ForeignOp' { arity: 1, op': op', op: fwd, op_bwd: unsafePa
    bwd :: OpBwd Unit
    bwd _ = error unimplemented
 
-type ArrayData a = Array2 (Val a)
-
 dims :: ForeignOp
 dims = mkExists $ ForeignOp' { arity: 1, op': op, op: fwd, op_bwd: unsafePartial bwd }
    where
@@ -108,12 +107,12 @@ dims = mkExists $ ForeignOp' { arity: 1, op': op, op: fwd, op_bwd: unsafePartial
       Constr <$> new (singleton α) <@> cPair <@> (v1 : v2 : Nil)
    op _ = throw "Matrix expected"
 
-   fwd :: OpFwd (Raw ArrayData)
+   fwd :: OpFwd (Array2 (Raw Val))
    fwd (Matrix α (MatrixRep (vss × (i × β1) × (j × β2))) : Nil) =
       pure $ (map erase <$> vss) × Constr α cPair (Int β1 i : Int β2 j : Nil)
    fwd _ = throw "Matrix expected"
 
-   bwd :: Partial => OpBwd (Raw ArrayData)
+   bwd :: Partial => OpBwd (Array2 (Raw Val))
    bwd (vss × Constr α c (Int β1 i : Int β2 j : Nil)) | c == cPair =
       Matrix α (MatrixRep (((<$>) botOf <$> vss) × (i × β1) × (j × β2))) : Nil
 
@@ -122,22 +121,52 @@ matrixLookup = mkExists $ ForeignOp' { arity: 2, op': op, op: fwd, op_bwd: bwd }
    where
    op :: OpGraph
    op (Matrix _ r : Constr _ c (Int _ i : Int _ j : Nil) : Nil)
-      | c == cPair = matrixGet i j r
+      | c == cPair = pure $ matrixGet i j r
    op _ = throw "Matrix and pair of integers expected"
 
-   fwd :: OpFwd (Raw ArrayData × (Int × Int) × (Int × Int))
+   fwd :: OpFwd (Raw MatrixRep × (Int × Int))
    fwd (Matrix _ r@(MatrixRep (vss × (i' × _) × (j' × _))) : Constr _ c (Int _ i : Int _ j : Nil) : Nil)
       | c == cPair = do
-           v <- matrixGet i j r
-           pure $ ((map erase <$> vss) × (i' × j') × (i × j)) × v
+           let v = matrixGet i j r
+           pure $ (MatrixRep ((map erase <$> vss) × ((i' × unit) × (j' × unit))) × (i × j)) × v
    fwd _ = throw "Matrix and pair of integers expected"
 
-   bwd :: OpBwd (Raw ArrayData × (Int × Int) × (Int × Int))
-   bwd ((vss × (i' × j') × (i × j)) × v) =
+   bwd :: OpBwd (Raw MatrixRep × (Int × Int))
+   bwd (((MatrixRep (vss × (i' × _) × (j' × _))) × (i × j)) × v) =
       Matrix bot (matrixUpdate i j (const v) (MatrixRep (((<$>) botOf <$> vss) × (i' × bot) × (j' × bot))))
          : Constr bot cPair (Int bot i : Int bot j : Nil)
          : Nil
 
+matrixMut :: ForeignOp
+matrixMut = mkExists $ ForeignOp' { arity: 3, op': op, op: fwd, op_bwd: bwd }
+   where
+   op :: OpGraph
+   op (Matrix α r : Constr _ c (Int _ i : Int _ j : Nil) : v : Nil)
+      | c == cPair = Matrix <$> new (singleton α) <@> (matrixUpdate i j (const v) r)
+   op _ = throw "Matrix, pair of ints, and new val expected"
+
+   fwd :: OpFwd (Raw MatrixRep × (Int × Int) × Raw Val)
+   fwd (Matrix α r@(MatrixRep (vss × (i' × _) × (j' × _))) : Constr _ c (Int _ i : Int _ j : Nil) : v : Nil)
+      | c == cPair =
+           let
+              oldV = matrixGet i j r
+              newM = matrixUpdate i j (const v) r
+           in
+              pure $ (MatrixRep ((map erase <$> vss) × ((i' × unit) × (j' × unit))) × (i × j) × (erase oldV)) × (Matrix α newM)
+
+   fwd _ = throw "Matrix, pair of ints, and new val expected"
+
+   bwd :: OpBwd (Raw MatrixRep × (Int × Int) × Raw Val)
+   bwd ((((MatrixRep (vss × (i' × _) × (j' × _))) × (i × j) × oldV) × (Matrix α r))) =
+      let
+         newV = matrixGet i j r
+      in
+         Matrix α (matrixUpdate i j (const (map (const bot) oldV)) (MatrixRep (((<$>) botOf <$> vss) × (i' × bot) × (j' × bot))))
+            : Constr bot cPair (Int bot i : Int bot j : Nil)
+            : newV
+            : Nil
+   bwd _ = error "absurd backwards!"
+
 dict_difference :: ForeignOp
 dict_difference = mkExists $ ForeignOp' { arity: 2, op': op, op: fwd, op_bwd: unsafePartial bwd }
    where

src/Val.purs

@@ -23,7 +23,7 @@ import Foreign.Object (keys) as O
 import Graph (Vertex(..))
 import Graph.GraphWriter (class MonadWithGraphAlloc)
 import Lattice (class BoundedJoinSemilattice, class BoundedLattice, class Expandable, class JoinSemilattice, class Neg, Raw, definedJoin, expand, maybeJoin, neg, (∨))
-import Util (type (×), Endo, error, orElse, throw, unsafeUpdateAt, (!), (×), (≜), (≞))
+import Util (type (×), Endo, definitely, error, orElse, throw, unsafeUpdateAt, (!), (×), (≜), (≞))
 import Util.Pretty (Doc, beside, text)
 
 data Val a
@@ -106,11 +106,10 @@ newtype MatrixRep a = MatrixRep (Array2 (Val a) × (Int × a) × (Int × a))
 
 type Array2 a = Array (Array a)
 
-matrixGet :: forall a m. MonadThrow Error m => Int -> Int -> MatrixRep a -> m (Val a)
-matrixGet i j (MatrixRep (vss × _ × _)) =
-   orElse "Index out of bounds" do
-      us <- vss !! (i - 1)
-      us !! (j - 1)
+matrixGet :: forall a. Int -> Int -> MatrixRep a -> Val a
+matrixGet i j (MatrixRep (vss × _ × _)) = definitely "index out of bounds!" $ do
+   us <- vss !! (i - 1)
+   us !! (j - 1)
 
 matrixUpdate :: forall a. Int -> Int -> Endo (Val a) -> Endo (MatrixRep a)
 matrixUpdate i j δv (MatrixRep (vss × h × w)) =

test/Main.purs

@@ -23,8 +23,20 @@ tests = concat
    , test_linking
    ]
 
--- test_scratchpad :: Array (String × Aff Unit)
--- test_scratchpad = []
+-- tests = test_scratchpad
+
+test_scratchpad :: Array (String × Aff Unit)
+test_scratchpad = second void <$> many
+   [ { file: "dtw/matrix-update"
+     , fwd_expect:
+          "FNum 0, Infty, Infty, Infty, Infty, Infty, Infty,\n\
+          \Infty, FNum 0, FNum 0, Infty, Infty, Infty, Infty,\n\
+          \Infty, FNum 0, FNum 0, FNum 0, Infty, Infty, Infty,\n\
+          \Infty, Infty, FNum 0, FNum 0, FNum 0, Infty, Infty,\n\
+          \Infty, Infty, Infty, FNum 0, FNum 0, FNum 0, Infty"
+     }
+   ]
+   1
 
 test_desugaring :: Array (String × Aff Unit)
 test_desugaring = second void <$> many desugar_cases 1

test/Spec/Specs.purs

@@ -40,6 +40,21 @@ misc_cases =
    , { file: "range", fwd_expect: "((0, 0) : ((0, 1) : ((1, 0) : ((1, 1) : []))))" }
    , { file: "records", fwd_expect: "{a : 2, b : 6, c : 7, d : (5 : []), e : 7}" }
    , { file: "reverse", fwd_expect: "(2 : (1 : []))" }
+   , { file: "dtw/next-indices", fwd_expect: "((1, 1) : ((1, 2) : ((1, 3) : ((2, 1) : ((2, 2) : ((2, 3) : ((2, 4) : ((3, 1) : ((3, 2) : ((3, 3) : ((3, 4) : ((3, 5) : ((4, 2) : ((4, 3) : ((4, 4) : ((4, 5) : ((4, 6) : ((5, 3) : ((5, 4) : ((5, 5) : ((5, 6) : ((5, 7) : []))))))))))))))))))))))" }
+   , { file: "dtw/cost-matrix"
+     , fwd_expect:
+          "FNum 0, Infty, Infty, Infty, Infty, Infty, Infty,\n\
+          \Infty, FNum 0, FNum 0, Infty, Infty, Infty, Infty,\n\
+          \Infty, FNum 0, FNum 0, FNum 0, Infty, Infty, Infty,\n\
+          \Infty, Infty, FNum 0, FNum 0, FNum 0, Infty, Infty,\n\
+          \Infty, Infty, Infty, FNum 0, FNum 0, FNum 0, Infty"
+     }
+   , { file: "dtw/matrix-update"
+     , fwd_expect:
+          "100, (1, 2), (1, 3),\n\
+          \(2, 1), (2, 2), (2, 3),\n\
+          \(3, 1), (3, 2), (3, 3)"
+     }
    ]
 
 desugar_cases :: Array TestSpec