diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,6 @@ +module Main (main) where + +import Distribution.Simple + +main :: IO () +main = defaultMain diff --git a/dtw.cabal b/dtw.cabal new file mode 100644 --- /dev/null +++ b/dtw.cabal @@ -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 diff --git a/src/Data/DTW.hs b/src/Data/DTW.hs new file mode 100644 --- /dev/null +++ b/src/Data/DTW.hs @@ -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 #-} diff --git a/test/MainTest.hs b/test/MainTest.hs new file mode 100644 --- /dev/null +++ b/test/MainTest.hs @@ -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-} + ]