packages feed

dtw (empty) → 0.9.2.0

raw patch · 5 files changed

+386/−0 lines, 5 filesdep +MemoTriedep +QuickCheckdep +basesetup-changed

Dependencies added: MemoTrie, QuickCheck, base, containers, test-framework, test-framework-quickcheck2, thyme, vector, vector-space

Files

+ LICENSE view
@@ -0,0 +1,20 @@+Copyright (c) 2014 Florian Hofmann++Permission is hereby granted, free of charge, to any person obtaining+a copy of this software and associated documentation files (the+"Software"), to deal in the Software without restriction, including+without limitation the rights to use, copy, modify, merge, publish,+distribute, sublicense, and/or sell copies of the Software, and to+permit persons to whom the Software is furnished to do so, subject to+the following conditions:++The above copyright notice and this permission notice shall be included+in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ Setup.hs view
@@ -0,0 +1,6 @@+module Main (main) where++import Distribution.Simple++main :: IO ()+main = defaultMain
+ dtw.cabal view
@@ -0,0 +1,61 @@+-- Initial dtw.cabal generated by cabal init.  For further +-- documentation, see http://haskell.org/cabal/users-guide/++name:                dtw+version:             0.9.2.0+synopsis:            (Fast) Dynamic Time Warping+description:         This package implements dynamic time warping as described+                     here http://en.wikipedia.org/w/index.php?title=Dynamic_time_warping&oldid=643501828+license:             MIT+license-file:        LICENSE+author:              Florian Hofmann+maintainer:          fho@f12n.de+-- copyright:           +category:            Algorithms, Data mining, Machine learning+build-type:          Simple+-- extra-source-files:  +cabal-version:       >=1.10++source-repository head+  type:                git+  location:            https://github.com/fhaust/dtw+++library+  build-depends:       base         >= 4.6 && < 4.9,+                       vector       >= 0.10 && < 0.11,+                       vector-space >= 0.10 && < 0.11,+                       containers   >= 0.5 && < 0.6,+                       MemoTrie     >= 0.6 && < 0.7+  exposed-modules:     Data.DTW+  ghc-options:         -Wall+  hs-source-dirs:      src+  default-language:    Haskell2010++test-suite maintest+  type:                exitcode-stdio-1.0+  hs-source-dirs:      src,test+  main-is:             MainTest.hs+  build-depends:       base,+                       vector,+                       vector-space,+                       thyme,+                       containers,+                       test-framework,+                       test-framework-quickcheck2,+                       QuickCheck <2.8,+                       MemoTrie+  ghc-options:         -O2 -Wall+  hs-source-dirs:      src+  default-language:    Haskell2010++--executable glossaer+--  main-is:             GlossAER.hs+--  -- other-modules:       +--  -- other-extensions:    +--  build-depends:       base >=4.7 && <4.8,+--                       aer,+--                       kdtree+--  ghc-options:         -O2 -Wall -fexcess-precision+--  hs-source-dirs:      src+--  default-language:    Haskell2010
+ src/Data/DTW.hs view
@@ -0,0 +1,221 @@+{-# LANGUAGE TypeFamilies #-}++-- | This module implements dynamic time warping as described here:+-- http://en.wikipedia.org/w/index.php?title=Dynamic_time_warping&oldid=643501828+--+-- Additionally 'fastDtw' is implemented as described in the paper:+-- "FastDTW: Toward Accurate Dynamic Time Warping in Linear Time and+-- Space" by Stan Salvador and Philip Chan.+--+-- Please note that 'fastDtw' is only an approximative solution. If you+-- need the optimal solution and can bear with the heavily increased demand+-- both in cpu and in memory you should use 'dtwMemo' or 'dtwMemoWindowed'.+--+--+-- == Example+--+-- >>> -- create two sample datasets+-- >>> let as = [ sin x | x <- [0,0.1..pi] ]+-- >>> let bs = [ sin (x+0.1) | x <- [0,0.1..pi] ]+-- >>> -- define a cost function between two datapoints+-- >>> let dist x y = abs (x-y)+-- >>> -- define a function that will half the size of a dataset (see below)+-- >>> let shrink xs = case xs of (a:b:cs) -> (a+b)/2 : shrink cs; a:[] -> [a]; [] -> []+-- >>> -- calculate the cost with fastDtw and dtwMemo for comparison+-- >>> cost $ fastDtw dist shrink 2 as bs :: Float+-- 0.19879311+-- >>> cost $ dtwMemo (\x y -> abs (x-y)) as bs :: Float+-- 0.19879311+--+-- == Some words on the shrink function+--+-- Care must be taken when choosing a shrink function. It's vital that the+-- resolution is halfed, this is not exactly a problem with the algorithm+-- but with the implementation. The lower resolution dataset should be an+-- honest representation of the higher resolution dataset. For starters+-- binning as in the example above should suffice.+-- +++module Data.DTW (dtwNaive, dtwMemo, fastDtw, DataSet(..), Result(..), Path, Index) where+++import qualified Data.Sequence as S++import qualified Data.Set as Set+import qualified Data.List as L+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as UV+import qualified Data.Vector.Storable as SV++import           Data.MemoTrie+import           Data.Function++++-- | a generic dataset is basically just an indexing function+-- | and an indicator of the dataset size+class DataSet dataset  where+    type Item dataset :: *+    ix  :: dataset -> Int -> Item dataset+    len :: dataset -> Int++-- some DataSet orphan instances+instance DataSet (S.Seq a) where+    type Item (S.Seq a) = a+    ix  = S.index+    len = S.length++instance DataSet [a] where+    type Item [a] = a+    ix  = (!!)+    len = length++instance DataSet (V.Vector a) where+    type Item (V.Vector a) = a+    ix  = V.unsafeIndex -- for speed?+    len = V.length++instance UV.Unbox a => DataSet (UV.Vector a) where+    type Item (UV.Vector a) = a+    ix  = UV.unsafeIndex -- for speed?+    len = UV.length++instance SV.Storable a => DataSet (SV.Vector a) where+    type Item (SV.Vector a) = a+    ix  = SV.unsafeIndex -- for speed?+    len = SV.length++-- common types++type Index  = (Int,Int)+type Path   = [Index]+type Window = Set.Set Index++data Result a = Result { cost :: a, path :: Path } deriving (Show,Read,Eq)+++-- | this is the naive implementation of dynamic time warping+-- no caching what so ever is taking place+-- this should not be used and is just used as a reference for the other+-- implementations++dtwNaive :: (Ord c, Fractional c, DataSet a, DataSet b)+         => (Item a -> Item b -> c) -> a -> b -> c+dtwNaive δ as bs = go (len as - 1) (len bs - 1)+    where go 0 0 = 0+          go _ 0 = 1/0+          go 0 _ = 1/0+          go x y = δ (ix as x) (ix bs y) + minimum [ go (x-1)  y+                                                   , go  x    (y-1)+                                                   , go (x-1) (y-1)+                                                   ]++-------------------------------------------------------------------------------------++-- | this is the "standard" implementation of dynamic time warping+-- O(N^2) is achieved by memoization of previous results+dtwMemo :: (Ord c, Fractional c, DataSet a, DataSet b)+        => (Item a -> Item b -> c) -> a -> b -> Result c+dtwMemo δ = dtwMemoWindowed δ (\_ _ -> True)++{-# INLINABLE dtwMemo #-}++-- | "standard" implementation of dynamic time warping with an additional+-- parameter that can be used to define a search window+dtwMemoWindowed :: (Ord c, Fractional c, DataSet a, DataSet b)+                => (Item a -> Item b -> c)+                -> (Int -> Int -> Bool)+                -> a+                -> b+                -> Result c+dtwMemoWindowed δ inWindow as bs = go (len as - 1) (len bs - 1)+    where -- wrap go' in a memoziation function so that each value+          -- is calculated only once+          go = memo2 go'+          -- handle special cases, origin cost is zero,+          -- border cost is infinity+          go' 0 0                      = Result 0 [(0,0)]+          go' 0 y                      = Result (1/0) [(0,y)]+          go' x 0                      = Result (1/0) [(x,0)]+          -- check that this index is not out of the search window+          go' x y | not (inWindow x y) = Result (1/0) [(x,y)]+          -- else calculate this value, note that this calls the+          -- memoized version of go recursivly+          go' x y                      = Result newCost newPath+            where minResult = L.minimumBy (compare `on` cost) [ go (x-1)  y+                                                              , go  x    (y-1)+                                                              , go (x-1) (y-1) ]+                  newPath   = (x,y) : path minResult+                  newCost   = δ (ix as x) (ix bs y) + cost minResult++{-# INLINABLE dtwMemoWindowed #-}++-------------------------------------------------------------------------------------++{--- | reduce a dataset to half its size by averaging neighbour values-}+{--- together-}+{-reduceByHalf :: Fractional a => Seq a -> Seq a-}+{-reduceByHalf (S.viewl -> x :< (S.viewl -> y :< xs))  = (x + y) / 2 <| reduceByHalf xs-}+{-reduceByHalf (S.viewl -> x :< (S.viewl -> S.EmptyL)) = S.singleton x-}+{-reduceByHalf _                                       = S.empty-}++-- | create a search window by projecting the path from a lower resolution+-- to the next level as defined in the fastdtw paper (figure 6), ie:+--+-- +-------+-------+        +---+---+---+---++-- |       |       |        |   |   | X | X |+-- |       |   X   |        +---+---+---+---++-- |       |       |        |   | X | X | X |+-- +-------+-------+   ->   +---+---+---+---++-- |       |       |        | X | X | X |   |+-- |   X   |       |        +---+---+---+---++-- |       |       |        | X | X |   |   |+-- +-------+-------+        +---+---+---+---++projectPath :: Path -> Window+projectPath p = Set.fromList $ concatMap expand $ concat $ zipWith project p (tail p)+  where project (a,b) (c,d) = [(2*a,2*b),((2*a+2*c) `quot` 2, (2*b+2*d) `quot` 2), (2*c,2*d)]+        expand  (a,b)       = [(a,b),(a+1,b),(a,b+1),(a+1,b+1)]++-- | expand the search window by a given radius+-- (compare fastdtw paper figure 6)+-- ie for a radius of 1:+--+-- +---+---+---+---+        +---+---+---+---++-- |   |   | X | X |        | o | o | X | X |+-- +---+---+---+---+        +---+---+---+---++-- |   | X | X | X |        | o | X | X | X |+-- +---+---+---+---+   ->   +---+---+---+---++-- | X | X | X |   |        | X | X | X | o |+-- +---+---+---+---+        +---+---+---+---++-- | X | X |   |   |        | X | X | o | o |+-- +---+---+---+---+        +---+---+---+---++--+-- the "o"s mark the expanded regions+expandWindow :: Int -> Window -> Window+expandWindow r = Set.fromList . concatMap (\(x,y) -> [ (x',y') | y' <- [y-r..y+r], x' <- [x-r..x+r] ]) . Set.toList++++-- | this is the "fast" implementation of dynamic time warping+-- as per the authors this methods calculates a good approximate+-- result in O(N), depending on the usecase the windowsize should be+-- tweaked++fastDtw :: (Ord c, Fractional c, DataSet a)+        => (Item a -> Item a -> c)+        -> (a -> a) -- ^ function that shrinks a dataset by a factor of two+        -> Int      -- ^ radius that the search window is expanded at each resolution level+        -> a        -- ^ first dataset+        -> a        -- ^ second dataset+        -> Result c -- ^ result+fastDtw δ shrink r as bs | len as <= minTSsize || len bs <= minTSsize = dtwMemo δ as bs+                         | otherwise = dtwMemoWindowed δ inWindow as bs+    where minTSsize    = r+2+          shrunkAS     = shrink as+          shrunkBS     = shrink bs+          lowResResult = fastDtw δ shrink r shrunkAS shrunkBS+          window       = expandWindow r $ projectPath (path lowResResult)+          inWindow x y = (x,y) `Set.member` window++{-# INLINABLE fastDtw #-}
+ test/MainTest.hs view
@@ -0,0 +1,78 @@++{-# LANGUAGE ViewPatterns #-}++module Main where++++-- module under test+import Data.DTW++import Test.Framework+import Test.Framework.Providers.QuickCheck2++import Test.QuickCheck++import           Data.Sequence (Seq(..), ViewL(..), (<|))+import qualified Data.Sequence as S++import Data.Functor+++newtype SmallNonEmptySeq a = SmallNonEmptySeq { getSmallNonEmpty :: [a] }+    deriving (Show, Eq)++instance Arbitrary a => Arbitrary (SmallNonEmptySeq a) where+    arbitrary = SmallNonEmptySeq <$> listOf arbitrary `suchThat` (\l -> length l > 2 && length l < 10)++newtype MediumNonEmptySeq a = MediumNonEmptySeq { getMediumNonEmpty :: [a] }+    deriving (Show, Eq)++instance Arbitrary a => Arbitrary (MediumNonEmptySeq a) where+    arbitrary = MediumNonEmptySeq <$> listOf arbitrary `suchThat` (\l -> length l > 100 && length l < 1000)+++dist :: Double -> Double -> Double+dist x y = abs (x-y)++-- | reduce a dataset to half its size by averaging neighbour values+-- together+reduceByHalf :: Fractional a => Seq a -> Seq a+reduceByHalf (S.viewl -> x :< (S.viewl -> y :< xs))  = (x + y) / 2 <| reduceByHalf xs+reduceByHalf (S.viewl -> x :< (S.viewl -> S.EmptyL)) = S.singleton x+reduceByHalf _                                       = S.empty++{-testDTWVSDTWNaive :: (SmallNonEmptySeq Double, SmallNonEmptySeq Double) -> Bool-}+{-testDTWVSDTWNaive (la,lb) = abs (dtwNaive dist sa sb - dtw dist sa sb) < 0.01-}+{-  where sa = S.fromList $ getSmallNonEmpty la-}+{-        sb = S.fromList $ getSmallNonEmpty lb-}++testDTWMemoVSDTWNaive :: (SmallNonEmptySeq Double, SmallNonEmptySeq Double) -> Bool+testDTWMemoVSDTWNaive (la,lb) = abs (dtwNaive dist sa sb - cost (dtwMemo dist sa sb)) < 0.01+  where sa = S.fromList $ getSmallNonEmpty la+        sb = S.fromList $ getSmallNonEmpty lb++testFastDTWvsDTWNaive :: (SmallNonEmptySeq Double, SmallNonEmptySeq Double) -> Bool+testFastDTWvsDTWNaive (la,lb) = abs (1 - (ca/l) / (cb/l)) < 0.1+  where sa = S.fromList $ getSmallNonEmpty la+        sb = S.fromList $ getSmallNonEmpty lb+        l  = fromIntegral $ S.length sa + S.length sb+        ca = dtwNaive dist sa sb+        cb = cost $ fastDtw dist reduceByHalf 2 sa sb++-- FIXME no real idea how to compare an optimal and an approximative+-- algorithm ... best bet below, but still failing tests+{-testFastDTWvsDTWMemo :: (MediumNonEmptySeq Double, MediumNonEmptySeq Double) -> Bool-}+{-testFastDTWvsDTWMemo (la,lb) = abs (1 - ((costA / matSize) / (costB / matSize))) < 0.1-}+{-  where sa = S.fromList $ getMediumNonEmpty la-}+{-        sb = S.fromList $ getMediumNonEmpty lb-}+{-        costA = cost (dtwMemo dist sa sb)-}+{-        costB = cost (fastDtw dist 10 sa sb)-}+{-        matSize = fromIntegral $ S.length sa * S.length sb-}++main :: IO ()+main = defaultMain +     [ testProperty "dtwMemo ≡ dtwNaive" testDTWMemoVSDTWNaive+     , testProperty "fastDtw ≅ dtwNaive" testFastDTWvsDTWNaive+     {-, testProperty "fastDtw == dtwMemo"  testFastDTWvsDTWMemo-}+     ]