Zapamiętywanie w Haskell?

Wszelkie wskazówki, jak skutecznie rozwiązać następującą funkcję w Haskellu dla dużych liczb (n > 108)

f(n) = max(n, f(n/2) + f(n/3) + f(n/4))

Widziałem przykłady zapamiętywania w Haskellu w celu rozwiązania liczb Fibonacciego, które obejmowały (leniwie) obliczanie wszystkich liczb Fibonacciego do wymaganego n. Ale w tym przypadku dla danego n wystarczy obliczyć bardzo niewiele wyników pośrednich.

Dzięki

Możemy to zrobić bardzo efektywnie, tworząc strukturę, którą możemy indeksować w czasie nieliniowym.

Ale najpierw,

{-# LANGUAGE BangPatterns #-}

import Data.Function (fix)

Zdefiniujmy f, ale niech używa „otwartej rekursji” zamiast wywoływać się bezpośrednio.

f :: (Int -> Int) -> Int -> Int
f mf 0 = 0
f mf n = max n $ mf (n `div` 2) +
                 mf (n `div` 3) +
                 mf (n `div` 4)

Możesz uzyskać niezapamiętać f, używającfix f

Pozwoli ci to sprawdzić, fczy to, co masz na myśli, dla małych wartości f, dzwoniąc, na przykład:fix f 123 = 144

Moglibyśmy to zapamiętać, definiując:

f_list :: [Int]
f_list = map (f faster_f) [0..]

faster_f :: Int -> Int
faster_f n = f_list !! n

To działa dość dobrze i zastępuje to, co miało zająć O (n ^ 3) czas, czymś, co zapamiętuje wyniki pośrednie.

Ale wciąż potrzeba czasu liniowego, aby indeksować, aby znaleźć zapamiętaną odpowiedź mf. Oznacza to, że wyniki takie jak:

*Main Data.List> faster_f 123801
248604

są do zaakceptowania, ale wynik nie skaluje się dużo lepiej. Możemy zrobić lepiej!

Najpierw zdefiniujmy nieskończone drzewo:

data Tree a = Tree (Tree a) a (Tree a)
instance Functor Tree where
    fmap f (Tree l m r) = Tree (fmap f l) (f m) (fmap f r)

Następnie zdefiniujemy sposób indeksowania do niego, abyśmy mogli zamiast tego znaleźć węzeł z indeksem nw czasie O (log n) :

index :: Tree a -> Int -> a
index (Tree _ m _) 0 = m
index (Tree l _ r) n = case (n - 1) `divMod` 2 of
    (q,0) -> index l q
    (q,1) -> index r q

... i możemy znaleźć drzewo pełne liczb naturalnych dla wygody, więc nie będziemy musieli bawić się tymi indeksami:

nats :: Tree Int
nats = go 0 1
    where
        go !n !s = Tree (go l s') n (go r s')
            where
                l = n + s
                r = l + s
                s' = s * 2

Ponieważ możemy indeksować, możesz po prostu przekonwertować drzewo na listę:

toList :: Tree a -> [a]
toList as = map (index as) [0..]

Możesz sprawdzić dotychczasową pracę, weryfikując, która toList natsdaje[0..]

Teraz,

f_tree :: Tree Int
f_tree = fmap (f fastest_f) nats

fastest_f :: Int -> Int
fastest_f = index f_tree

działa tak samo, jak w przypadku powyższej listy, ale zamiast potrzebować czasu liniowego na znalezienie każdego węzła, można go ścigać w czasie logarytmicznym.

Wynik jest znacznie szybszy:

*Main> fastest_f 12380192300
67652175206

*Main> fastest_f 12793129379123
120695231674999

W rzeczywistości jest to o wiele szybsze, że możesz przejść przez powyższe i zastąpić Intje Integerpowyższymi i niemal natychmiast uzyskać absurdalnie duże odpowiedzi

*Main> fastest_f' 1230891823091823018203123
93721573993600178112200489

*Main> fastest_f' 12308918230918230182031231231293810923
11097012733777002208302545289166620866358

Odpowiedź Edwarda jest tak wspaniałą perełką, że ją zduplikowałem i dostarczyłem implementacje memoListi memoTreekombinatory, które zapamiętują funkcję w formie otwarto-rekurencyjnej.

{-# LANGUAGE BangPatterns #-}

import Data.Function (fix)

f :: (Integer -> Integer) -> Integer -> Integer
f mf 0 = 0
f mf n = max n $ mf (div n 2) +
                 mf (div n 3) +
                 mf (div n 4)


-- Memoizing using a list

-- The memoizing functionality depends on this being in eta reduced form!
memoList :: ((Integer -> Integer) -> Integer -> Integer) -> Integer -> Integer
memoList f = memoList_f
  where memoList_f = (memo !!) . fromInteger
        memo = map (f memoList_f) [0..]

faster_f :: Integer -> Integer
faster_f = memoList f


-- Memoizing using a tree

data Tree a = Tree (Tree a) a (Tree a)
instance Functor Tree where
    fmap f (Tree l m r) = Tree (fmap f l) (f m) (fmap f r)

index :: Tree a -> Integer -> a
index (Tree _ m _) 0 = m
index (Tree l _ r) n = case (n - 1) `divMod` 2 of
    (q,0) -> index l q
    (q,1) -> index r q

nats :: Tree Integer
nats = go 0 1
    where
        go !n !s = Tree (go l s') n (go r s')
            where
                l = n + s
                r = l + s
                s' = s * 2

toList :: Tree a -> [a]
toList as = map (index as) [0..]

-- The memoizing functionality depends on this being in eta reduced form!
memoTree :: ((Integer -> Integer) -> Integer -> Integer) -> Integer -> Integer
memoTree f = memoTree_f
  where memoTree_f = index memo
        memo = fmap (f memoTree_f) nats

fastest_f :: Integer -> Integer
fastest_f = memoTree f

Nie jest to najbardziej efektywny sposób, ale zapamiętuje:

f = 0 : [ g n | n <- [1..] ]
    where g n = max n $ f!!(n `div` 2) + f!!(n `div` 3) + f!!(n `div` 4)

przy żądaniu f !! 144sprawdza się, czy f !! 143istnieje, ale jego dokładna wartość nie jest obliczana. Nadal jest to nieznany wynik obliczeń. Jedynymi dokładnymi obliczonymi wartościami są te potrzebne.

Tak więc początkowo, o ile zostało obliczone, program nic nie wie.

f = .... 

Kiedy wysyłamy żądanie f !! 12, zaczyna ono dopasowywać wzorce:

f = 0 : g 1 : g 2 : g 3 : g 4 : g 5 : g 6 : g 7 : g 8 : g 9 : g 10 : g 11 : g 12 : ...

Teraz zaczyna obliczać

f !! 12 = g 12 = max 12 $ f!!6 + f!!4 + f!!3

To rekurencyjnie tworzy kolejne żądanie na f, więc obliczamy

f !! 6 = g 6 = max 6 $ f !! 3 + f !! 2 + f !! 1
f !! 3 = g 3 = max 3 $ f !! 1 + f !! 1 + f !! 0
f !! 1 = g 1 = max 1 $ f !! 0 + f !! 0 + f !! 0
f !! 0 = 0

Teraz możemy trochę cofnąć

f !! 1 = g 1 = max 1 $ 0 + 0 + 0 = 1

Co oznacza, że ​​program teraz wie:

f = 0 : 1 : g 2 : g 3 : g 4 : g 5 : g 6 : g 7 : g 8 : g 9 : g 10 : g 11 : g 12 : ...

Nadal sączy się:

f !! 3 = g 3 = max 3 $ 1 + 1 + 0 = 3

Co oznacza, że ​​program teraz wie:

f = 0 : 1 : g 2 : 3 : g 4 : g 5 : g 6 : g 7 : g 8 : g 9 : g 10 : g 11 : g 12 : ...

Teraz kontynuujemy obliczenia f!!6:

f !! 6 = g 6 = max 6 $ 3 + f !! 2 + 1
f !! 2 = g 2 = max 2 $ f !! 1 + f !! 0 + f !! 0 = max 2 $ 1 + 0 + 0 = 2
f !! 6 = g 6 = max 6 $ 3 + 2 + 1 = 6

Co oznacza, że ​​program teraz wie:

f = 0 : 1 : 2 : 3 : g 4 : g 5 : 6 : g 7 : g 8 : g 9 : g 10 : g 11 : g 12 : ...

Teraz kontynuujemy obliczenia f!!12:

f !! 12 = g 12 = max 12 $ 6 + f!!4 + 3
f !! 4 = g 4 = max 4 $ f !! 2 + f !! 1 + f !! 1 = max 4 $ 2 + 1 + 1 = 4
f !! 12 = g 12 = max 12 $ 6 + 4 + 3 = 13

Co oznacza, że ​​program teraz wie:

f = 0 : 1 : 2 : 3 : 4 : g 5 : 6 : g 7 : g 8 : g 9 : g 10 : g 11 : 13 : ...

Więc obliczenia są wykonywane dość leniwie. Program wie, że f !! 8istnieje jakaś wartość for , że jest równa g 8, ale nie ma pojęcia, co nią g 8jest.

To jest dodatek do doskonałej odpowiedzi Edwarda Kmetta.

Kiedy wypróbowałem jego kod, definicje natsi indexwydawały się dość tajemnicze, więc piszę alternatywną wersję, która była dla mnie łatwiejsza do zrozumienia.

Definiuję indexi natsw kategoriach index'i nats'.

index' t njest zdefiniowany w zakresie [1..]. (Przypomnij sobie, że index tjest to zdefiniowane w zakresie [0..].) To działa przeszukuje drzewo traktując je njako ciąg bitów i odczytując je w odwrotnej kolejności. Jeśli bit jest 1, bierze prawą gałąź. Jeśli bit jest 0, bierze lewą gałąź. Zatrzymuje się, gdy osiągnie ostatni bit (który musi być a 1).

index' (Tree l m r) 1 = m
index' (Tree l m r) n = case n `divMod` 2 of
                          (n', 0) -> index' l n'
                          (n', 1) -> index' r n'

Tak jak natsjest zdefiniowane dla, indextak index nats n == njest zawsze prawdziwe, nats'jest zdefiniowane dla index'.

nats' = Tree l 1 r
  where
    l = fmap (\n -> n*2)     nats'
    r = fmap (\n -> n*2 + 1) nats'
    nats' = Tree l 1 r

Teraz natsi indexsą po prostu nats'i index'ale z wartościami przesuniętymi o 1:

index t n = index' t (n+1)
nats = fmap (\n -> n-1) nats'

Jak stwierdzono w odpowiedzi Edwarda Kmetta, aby przyspieszyć działanie, trzeba buforować kosztowne obliczenia i mieć do nich szybki dostęp.

Aby funkcja nie była monadyczna, rozwiązanie polegające na budowaniu nieskończonego leniwego drzewa z odpowiednim sposobem jego indeksowania (jak pokazano w poprzednich postach) spełnia ten cel. Jeśli zrezygnujesz z niemonadycznego charakteru funkcji, możesz użyć standardowych kontenerów asocjacyjnych dostępnych w Haskell w połączeniu z monadami „stanowymi” (takimi jak State lub ST).

Chociaż główną wadą jest to, że otrzymujesz funkcję niemonadyczną, nie musisz już samodzielnie indeksować struktury i możesz po prostu użyć standardowych implementacji kontenerów asocjacyjnych.

Aby to zrobić, musisz najpierw ponownie napisać swoją funkcję, aby akceptowała jakąkolwiek monadę:

fm :: (Integral a, Monad m) => (a -> m a) -> a -> m a
fm _    0 = return 0
fm recf n = do
   recs <- mapM recf $ div n <$> [2, 3, 4]
   return $ max n (sum recs)

W swoich testach nadal możesz zdefiniować funkcję, która nie wykonuje zapamiętywania za pomocą Data.Function.fix, chociaż jest nieco bardziej szczegółowa:

noMemoF :: (Integral n) => n -> n
noMemoF = runIdentity . fix fm

Następnie możesz użyć State monad w połączeniu z Data.Map, aby przyspieszyć działanie:

import qualified Data.Map.Strict as MS

withMemoStMap :: (Integral n) => n -> n
withMemoStMap n = evalState (fm recF n) MS.empty
   where
      recF i = do
         v <- MS.lookup i <$> get
         case v of
            Just v' -> return v' 
            Nothing -> do
               v' <- fm recF i
               modify $ MS.insert i v'
               return v'

Po drobnych zmianach można zamiast tego dostosować kod do pracy z Data.HashMap:

import qualified Data.HashMap.Strict as HMS

withMemoStHMap :: (Integral n, Hashable n) => n -> n
withMemoStHMap n = evalState (fm recF n) HMS.empty
   where
      recF i = do
         v <- HMS.lookup i <$> get
         case v of
            Just v' -> return v' 
            Nothing -> do
               v' <- fm recF i
               modify $ HMS.insert i v'
               return v'

Zamiast trwałych struktur danych możesz także wypróbować zmienne struktury danych (takie jak Data.HashTable) w połączeniu z monadą ST:

import qualified Data.HashTable.ST.Linear as MHM

withMemoMutMap :: (Integral n, Hashable n) => n -> n
withMemoMutMap n = runST $
   do ht <- MHM.new
      recF ht n
   where
      recF ht i = do
         k <- MHM.lookup ht i
         case k of
            Just k' -> return k'
            Nothing -> do 
               k' <- fm (recF ht) i
               MHM.insert ht i k'
               return k'

W porównaniu z implementacją bez zapamiętywania, każda z tych implementacji pozwala, przy dużych nakładach, na uzyskanie wyników w mikro-sekundach, zamiast czekać kilka sekund.

Używając Criterion jako benchmarku, mogłem zauważyć, że implementacja z Data.HashMap faktycznie działała nieco lepiej (około 20%) niż ta z Data.Map i Data.HashTable, dla których czasy były bardzo podobne.

Wyniki testu porównawczego okazały się nieco zaskakujące. Moje początkowe wrażenie było takie, że HashTable przewyższy implementację HashMap, ponieważ jest zmienna. W tej ostatniej implementacji może być ukryty jakiś defekt wydajności.

Kilka lat później przyjrzałem się temu i zdałem sobie sprawę, że istnieje prosty sposób na zapamiętanie tego w czasie liniowym przy użyciu zipWithfunkcji pomocniczej:

dilate :: Int -> [x] -> [x]
dilate n xs = replicate n =<< xs

dilatema tę przydatną właściwość dilate n xs !! i == xs !! div i n.

Więc zakładając, że mamy f (0), to upraszcza obliczenia do

fs = f0 : zipWith max [1..] (tail $ fs#/2 .+. fs#/3 .+. fs#/4)
  where (.+.) = zipWith (+)
        infixl 6 .+.
        (#/) = flip dilate
        infixl 7 #/

Wygląda bardzo podobnie do naszego oryginalnego opisu problemu i daje liniowe rozwiązanie ( sum $ take n fszajmie O (n)).

Kolejny dodatek do odpowiedzi Edwarda Kmetta: samodzielny przykład:

data NatTrie v = NatTrie (NatTrie v) v (NatTrie v)

memo1 arg_to_index index_to_arg f = (\n -> index nats (arg_to_index n))
  where nats = go 0 1
        go i s = NatTrie (go (i+s) s') (f (index_to_arg i)) (go (i+s') s')
          where s' = 2*s
        index (NatTrie l v r) i
          | i <  0    = f (index_to_arg i)
          | i == 0    = v
          | otherwise = case (i-1) `divMod` 2 of
             (i',0) -> index l i'
             (i',1) -> index r i'

memoNat = memo1 id id 

Użyj go w następujący sposób, aby zapamiętać funkcję z pojedynczą liczbą całkowitą arg (np. Fibonacci):

fib = memoNat f
  where f 0 = 0
        f 1 = 1
        f n = fib (n-1) + fib (n-2)

Buforowane będą tylko wartości argumentów nieujemnych.

Aby również buforować wartości argumentów ujemnych, użyj memoInt, zdefiniowane w następujący sposób:

memoInt = memo1 arg_to_index index_to_arg
  where arg_to_index n
         | n < 0     = -2*n
         | otherwise =  2*n + 1
        index_to_arg i = case i `divMod` 2 of
           (n,0) -> -n
           (n,1) ->  n

Aby buforować wartości dla funkcji z dwoma argumentami całkowitymi, użyj memoIntInt, zdefiniowanych w następujący sposób:

memoIntInt f = memoInt (\n -> memoInt (f n))

Rozwiązanie bez indeksowania, a nie oparte o Edward KMETT.

Oddzielam wspólne poddrzewa do wspólnego rodzica ( f(n/4)jest dzielone między f(n/2)i f(n/4)i f(n/6)jest dzielone między f(2)i f(3)). Zapisując je jako pojedynczą zmienną w rodzicu, obliczenie poddrzewa jest wykonywane raz.

data Tree a =
  Node {datum :: a, child2 :: Tree a, child3 :: Tree a}

f :: Int -> Int
f n = datum root
  where root = f' n Nothing Nothing


-- Pass in the arg
  -- and this node's lifted children (if any).
f' :: Integral a => a -> Maybe (Tree a) -> Maybe (Tree a)-> a
f' 0 _ _ = leaf
    where leaf = Node 0 leaf leaf
f' n m2 m3 = Node d c2 c3
  where
    d = if n < 12 then n
            else max n (d2 + d3 + d4)
    [n2,n3,n4,n6] = map (n `div`) [2,3,4,6]
    [d2,d3,d4,d6] = map datum [c2,c3,c4,c6]
    c2 = case m2 of    -- Check for a passed-in subtree before recursing.
      Just c2' -> c2'
      Nothing -> f' n2 Nothing (Just c6)
    c3 = case m3 of
      Just c3' -> c3'
      Nothing -> f' n3 (Just c6) Nothing
    c4 = child2 c2
    c6 = f' n6 Nothing Nothing

    main =
      print (f 123801)
      -- Should print 248604.

Kod nie jest łatwo rozszerzany do ogólnej funkcji zapamiętywania (przynajmniej nie wiedziałbym, jak to zrobić) i naprawdę musisz pomyśleć, w jaki sposób podproblemy się pokrywają, ale strategia powinna działać dla ogólnych parametrów wielu niecałkowitych . (Wymyśliłem to dla dwóch parametrów ciągów.)

Notatka jest odrzucana po każdym obliczeniu. (Ponownie myślałem o dwóch parametrach ciągów.)

Nie wiem, czy to jest bardziej wydajne niż inne odpowiedzi. Technicznie każde wyszukiwanie składa się z jednego lub dwóch kroków („Spójrz na swoje dziecko lub dziecko”), ale może wymagać dużo dodatkowej pamięci.

Edycja: to rozwiązanie nie jest jeszcze poprawne. Udostępnianie jest niekompletne.

Edycja: Powinno teraz poprawnie udostępniać poddzieci, ale zdałem sobie sprawę, że ten problem ma wiele nietrywialnych współdzielenia: n/2/2/2i n/3/3może być taki sam. Problem nie pasuje do mojej strategii.