-- | K-means algorithm
kmeans' :: G.Vector v a
        => MU.Matrix Double         -- ^ Initial set of k centroids
        -> Int                      -- ^ Max inter
        -> v a                      -- ^ Input data
        -> (a -> U.Vector Double)   -- ^ Feature extraction function
        -> (U.Vector Int, MU.Matrix Double, Double)
kmeans' initial maxiter dat fn
    | U.length (fn $ G.head dat) /= d = error "Dimension mismatched."
    | otherwise = (member, centers, U.sum $ U.imap ( \i x -> sqrt $ sumSquares
        (fn $ dat G.! i) (centers `MU.takeRow` x) ) member )
  where
    (member, centers) = loop 0 initial U.empty
    loop !iter means membership
        | iter >= maxiter || membership' == membership = (membership, means)
        | otherwise = loop (iter+1) (update membership') membership'
      where
        membership' = assign means

    -- Assignment step
    assign means = U.generate n $ \i ->
        let x = fn $ G.unsafeIndex dat i
            f (!min', !j') j = let d = sumSquares x $ means `unsafeTakeRow` j
                               in if d < min' then (d, j) else (min', j')
        in snd $ foldl' f (1/0, -1) [0..k-1]

    -- Update step
    update membership = MU.create $ do
        m <- MM.replicate (k,d) 0.0
        count <- UM.replicate k (0 :: Int)
        forM_ [0..n-1] $ \i -> do
            let x = membership `U.unsafeIndex` i
                vec = fn $ dat `G.unsafeIndex` i
            UM.unsafeModify count (+1) x
            forM_ [0..d-1] $ \j ->
                MM.unsafeRead m (x,j) >>=
                    MM.unsafeWrite m (x,j) . (+ (vec `U.unsafeIndex` j))
        -- normalize
        forM_ [0..k-1] $ \i -> do
            c <- UM.unsafeRead count i
            forM_ [0..d-1] $ \j ->
                MM.unsafeRead m (i,j) >>= MM.unsafeWrite m (i,j) . (/fromIntegral c)
        return m
    n = G.length dat
    k = MU.rows initial
    d = MU.cols initial

kmeans' :: G.Vector v a
        => MU.Matrix Double         -- ^ Initial set of k centroids
        -> Int                      -- ^ Max inter
        -> v a                      -- ^ Input data
        -> (a -> U.Vector Double)   -- ^ Feature extraction function
        -> (U.Vector Int, MU.Matrix Double, Double)
kmeans' initial maxiter dat fn
    | U.length (fn $ G.head dat) /= d = error "Dimension mismatched."
    | otherwise = (member, centers, U.sum $ U.imap ( \i x -> sqrt $ sumSquares
        (fn $ dat G.! i) (centers `MU.takeRow` x) ) member )
  where
    (member, centers) = loop 0 initial U.empty
    loop !iter means membership
        | iter >= maxiter || membership' == membership = (membership, means)
        | otherwise = loop (iter+1) (update membership') membership'
      where
        membership' = assign means

    -- Assignment step
    assign means = U.generate n $ \i ->
        let x = fn $ G.unsafeIndex dat i
            f (!min', !j') j = let d = sumSquares x $ means `unsafeTakeRow` j
                               in if d < min' then (d, j) else (min', j')
        in snd $ foldl' f (1/0, -1) [0..k-1]

    -- Update step
    update membership = MU.create $ do
        m <- MM.replicate (k,d) 0.0
        count <- UM.replicate k (0 :: Int)
        forM_ [0..n-1] $ \i -> do
            let x = membership `U.unsafeIndex` i
                vec = fn $ dat `G.unsafeIndex` i
            UM.unsafeModify count (+1) x
            forM_ [0..d-1] $ \j ->
                MM.unsafeRead m (x,j) >>=
                    MM.unsafeWrite m (x,j) . (+ (vec `U.unsafeIndex` j))
        -- normalize
        forM_ [0..k-1] $ \i -> do
            c <- UM.unsafeRead count i
            forM_ [0..d-1] $ \j ->
                MM.unsafeRead m (i,j) >>= MM.unsafeWrite m (i,j) . (/fromIntegral c)
        return m
    n = G.length dat
    k = MU.rows initial
    d = MU.cols initial

kmeans' initial maxiter dat fn
    | U.length (fn $ G.head dat) /= d = error "Dimension mismatched."
    | otherwise = (member, centers, U.sum $ U.imap ( \i x -> sqrt $ sumSquares
        (fn $ dat G.! i) (centers `MU.takeRow` x) ) member )
  where
    (member, centers) = loop 0 initial U.empty
    loop !iter means membership
        | iter >= maxiter || membership' == membership = (membership, means)
        | otherwise = loop (iter+1) (update membership') membership'
      where
        membership' = assign means

    -- Assignment step
    assign means = U.generate n $ \i ->
        let x = fn $ G.unsafeIndex dat i
            f (!min', !j') j = let d = sumSquares x $ means `unsafeTakeRow` j
                               in if d < min' then (d, j) else (min', j')
        in snd $ foldl' f (1/0, -1) [0..k-1]

    -- Update step
    update membership = MU.create $ do
        m <- MM.replicate (k,d) 0.0
        count <- UM.replicate k (0 :: Int)
        forM_ [0..n-1] $ \i -> do
            let x = membership `U.unsafeIndex` i
                vec = fn $ dat `G.unsafeIndex` i
            UM.unsafeModify count (+1) x
            forM_ [0..d-1] $ \j ->
                MM.unsafeRead m (x,j) >>=
                    MM.unsafeWrite m (x,j) . (+ (vec `U.unsafeIndex` j))
        -- normalize
        forM_ [0..k-1] $ \i -> do
            c <- UM.unsafeRead count i
            forM_ [0..d-1] $ \j ->
                MM.unsafeRead m (i,j) >>= MM.unsafeWrite m (i,j) . (/fromIntegral c)
        return m
    n = G.length dat
    k = MU.rows initial
    d = MU.cols initial
