В прошлой главе мы узнали, что такое ленивые вычисления. В этой главе мы посмотрим чем они хороши. С ними можно делать невозможные вещи. Обращаться к ещё не вычисленным значениям, работать с бесконечными данными.
Мы пишем программу, чтобы решить какую-нибудь сложную задачу. Часто так бывает, что сложная задача оказывается сложной до тех пор пока её не удаётся разбить на отдельные независимые подзадачи. Мы решаем задачи по-меньше, потом собираем из них решения, из этих решений собираем другие решения и вот уже готова программа. Но мы решаем задачу не на листочке, нам необходимо объяснить её компьютеру. И тот язык, на котором мы пишем программу, оказывает сильное влияние на то как мы будем решать задачу. Мы не можем разбить программу на независимые подзадачи, если в том языке на котором мы собираемся объяснять задачу компьютеру нет средств для того, чтобы собрать эти решения вместе.
Об этом говорит Джон Хьюз (John Huges) в статье “Why functional programming matters”. Он приводит такую метафору. Если мы делаем стул и у нас нет хорошего клея. Единственное что нам остаётся это вырезать из дерева стул целиком. Это невероятно трудная задача. Гораздо проще сделать отдельные части и потом собрать вместе. Функциональные языки программирования предоставляют два новых вида “клея”. Это функции высшего порядка и ленивые вычисления. В статье можно найти много примеров. Некоторые из них мы рассмотрим в этой главе.
С функциями высших порядков мы уже знакомы, они позволяют склеивать небольшие решения. С их помощью мы можем параметризовать функцию другой функцией (поведением). Они дают нам возможность выделять сложные закономерности и собирать их в функции. Ленивые вычисления же предназначены для склеивания больших программ. Они синхронизируют выполнение подзадач, избавляя нас от необходимости выполнять это вручную.
Эта идея разбиения программы на независимые части приводит нас к понятию модульности. Когда мы решаем задачу мы пытаемся разложить её на простейшие составляющие. При этом часто оказывается, что эти составляющие применимы не только для нашей задачи, но и для многих других. Мы получаем целый букет решений, там где искали одно.
Рассмотрим несколько численных методов. Все эти методы построены на понятии сходимости. У нас есть последовательность решений и она сходится к одному решению, но мы не знаем когда. Мы только знаем, что промежуточные решения будут всё ближе и ближе к итоговому.
Поскольку у нас ленивый язык мы сначала построим все возможные решения, а затем выберем итоговое. Так же как мы делали это в прошлой главе, когда искали корни уравнения методом неподвижной точки. Эти примеры взяты из статьи “Why functional programming matters” Джона Хьюза.
Найдём производную функции в точке. Посмотрим на математическое определение производной:
$$f^\prime(x)\ =\lim_{h\to0}\frac{f(x+h)-f(x)}{h}$$
Производная это предел последовательности таких отношений, при $h$ стремящемся к нулю. Если предел сходится, то производная определена. Для того чтобы решить эту задачу мы начнём с небольшого значения $h$ и будем постепенно уменьшать его, вычисляя промежуточные значения производной. Как только они перестанут сильно изменяться мы будем считать, что мы нашли предел последовательности
Этот процесс напоминает то, что мы делали при поиске корня уравнения методом неподвижной точки. Мы можем взять из того решения функцию определения сходимости последовательности:
converge :: (Ord a, Num a) => a -> [a] -> a
converge eps (a:b:xs)
| abs (a - b) <= eps = a
| otherwise = converge eps (b:xs)
Теперь осталось только создать последовательность значений производных. Напишем функцию, которая вычисляет промежуточные решения:
easydiff :: Fractional a => (a -> a) -> a -> a -> a easydiff f x h = (f (x + h) - f x) / h
Мы возьмём начальное значение шага и будем последовательно уменьшать его вдвое:
halves = iterate (/2)
Соберём все части вместе:
diff :: (Ord a, Fractional a) => a -> a -> (a -> a) -> a -> a
diff h0 eps f x = converge eps $ map (easydiff f x) $ iterate (/2) h0
where easydiff f x h = (f (x + h) - f x) / h
Сохраним эти определения в отдельном модуле и найдём производную какой-нибудь функции. Протестируем решение на экспоненте. Известно, что производная экспоненты равна самой себе:
*Numeric> let exp' = diff 1 1e-5 exp *Numeric> let test x = abs $ exp x - exp' x *Numeric> test 2 1.4093421286887065e-5 *Numeric> test 5 1.767240203776055e-5
Теперь давайте поинтегрируем функции одного аргумента. Интеграл это площадь кривой под графиком функции. Если бы кривая была прямой, то мы могли бы вычислить интеграл по формуле трапеций:
easyintegrate :: Fractional a => (a -> a) -> a -> a -> a easyintegrate f a b = (f a + f b) * (b - a) / 2
Но мы хотим интегрировать не только прямые линии. Мы представим, что функция является ломаной прямой линией. Мы посчитаем интеграл на каждом из участков и сложим ответы. При этом чем ближе точки друг к другу, тем точнее можно представить функцию в виде ломаной прямой линии, тем точнее будет значение интеграла.
Проблема в том, что мы не знаем заранее насколько близки должны быть точки друг к другу. Это зависит от функции, которую мы хотим проинтегрировать. Но мы можем построить последовательность решений. На каждом шаге мы будем приближать функцию ломаной прямой, и на каждом шаге число изломов будет расти вдвое. Как только решение перестанет меняться мы вернём ответ.
Построим последовательность решений:
integrate :: Fractional a => (a -> a) -> a -> a -> [a]
integrate f a b = easyintegrate f a b :
zipWith (+) (integrate a mid) (integrate mid b)
where mid = (a + b)/2
Первое решение является площадью под прямой, которая соединяет концы отрезка. Потом мы делим отрезок пополам, строим последовательность приближений и складываем частичные суммы с помощью функции zipWith.
Эта версия функции хоть и наглядная, но не эффективная. Функция f вычисляется заново при каждом рекурсивном вызове. Было бы хорошо вычислять её только для новых значений. Для этого мы будем передавать значения с предыдущего шага:
integrate :: Fractional a => (a -> a) -> a -> a -> [a]
integrate f a b = integ f a b (f a) (f b)
where integ f a b fa fb = (fa+fb)*(b-a)/2 :
zipWith (+) (integ f a m fa fm)
(integ f m b fm fb)
where m = (a + b)/2
fm = f m
В этой версии мы вычисляем значения в функции f лишь один раз для каждой точки. Запишем итоговое решение:
int :: (Ord a, Fractional a) => a -> (a -> a) -> a -> a -> a int eps f a b = converge eps $ integrate f a b
Мы опять воспользовались функцией converge, нам не нужно было её переписывать. Проверим решение. Для проверки также воспользуемся экспонентой. В прошлой главе мы узнали, что
$$e^x = 1 + \int_0^x e^t dt$$
Посмотрим, так ли это для нашего алгоритма:
*Numeric> let exp' = int 1e-5 exp 0 *Numeric> let test x = abs $ exp x - 1 - exp' x *Numeric> test 2 8.124102876649886e-6 *Numeric> test 5 4.576306736225888e-6 *Numeric> test 10 1.0683757864171639e-5
Алгоритм работает. В статье ещё рассмотрены методы повышения точности этих алгоритмов. Что интересно для улучшения точности не надо менять существующий код. Функция принимает последовательность промежуточных решений и преобразует её.
Напишем модуль для вычисления степенных рядов. Этот пример взят из статьи Дугласа МакИлроя (Douglas McIlroy) “Power Series, Power Serious”. Степенной ряд представляет собой функцию, которая определяется списком коэффициентов:
$$F(x) = f_0 + f_1 x + f_2 x^2 + f_3 x^3 + f_4 x^4 + ...$$
Степенной ряд содержит бесконечное число слагаемых. Для вычисления нам потребуются функции сложения и умножения. Ряд $F(x)$ можно записать и по-другому:
| $F(x)$ | $=$ | $F_0 (x)$ |
| $=$ | $f_0 + x F_1 (x)$ | |
| $=$ | $f_0 + x (f_1 + x F_2 (x))$ |
Это определение очень похоже на определение списка. Ряд есть коэффициент $f_0$ и другой ряд $F_1(x)$ умноженный на x. Поэтому для представления рядов мы выберем конструкцию похожую на список:
data Ps a = a :+: Ps a
deriving (Show, Eq)
Но в нашем случае списки бесконечны, поэтому у нас лишь один конструктор. Далее мы будем писать просто $f + x F_1$, без скобок для аргумента.
Определим вспомогательные функции для создания рядов:
p0 :: Num a => a -> Ps a p0 x = x :+: p0 0 ps :: Num a => [a] -> Ps a ps [] = p0 0 ps (a:as) = a :+: ps as
Обратите внимание на то, как мы дописываем бесконечный хвост нулей в конец ряда. Теперь давайте определим функцию вычисления ряда. Мы будем вычислять лишь конечное число степеней.
eval :: Num a => Int -> Ps a -> a -> a eval 0 _ _ = 0 eval n (a :+: p) x = a + x * eval (n-1) p x
В первом случае мы хотим вычислить ноль степеней ряда, поэтому мы возвращаем ноль, а во втором случае значение ряда $a + x P$ складывается из числа $a$ и значения ряда $P$ умноженного на заданное значение.
В результате сложения и умножения рядов также получается ряд. Также мы можем создать ряд из числа. Эти операции говорят о том, что мы можем сделать степенной ряд экземпляром класса Num.
Рекурсивное представление ряда $f + x F$ позволяет нам очень кратко выражать операции, которые мы хотим определить. Теперь у нас нет бесконечного набора коэффициентов, у нас всего лишь одно число и ещё один ряд. Операции существенно упрощаются. Так сложение двух бесконечных рядов имеет вид:
$$F + G = (f + x F_1) + (g + x G_1) = (f+g) + x (F_1 + G_1)$$
Переведём на Haskell:
(f :+: fs) + (g :+: gs) = (f + g) :+: (fs + gs)
Умножим два ряда:
$$F*G = (f + x F_1) * (g + x G_1) = f g + x (f G_1 + F_1 * G)$$
Переведём:
(.*) :: Num a => a -> Ps a -> Ps a k .* (f :+: fs) = (k * f) :+: (k .* fs) (f :+: fs) * (g :+: gs) = (f * g) :+: (f .* gs + fs * (g :+: gs))
Дополнительная операция (.*) выполняет умножение всех коэффициентов ряда на число.
Соберём определения для методов класса Num вместе:
instance Num a => Num (Ps a) where
(f :+: fs) + (g :+: gs) = (f + g) :+: (fs + gs)
(f :+: fs) * (g :+: gs) = (f * g) :+: (f .* gs + fs * (g :+: gs))
negate (f :+: fs) = negate f :+: negate fs
fromInteger n = p0 (fromInteger n)
(.*) :: Num a => a -> Ps a -> Ps a
k .* (f :+: fs) = (k * f) :+: (k .* fs)
Методы abs и signum не определены для рядов. Обратите внимание на то, как рекурсивное определение рядов приводит к рекурсивным определениям функций для рядов. Этот приём очень характерен для Haskell. Поскольку наш ряд это число и ещё один ряд за счёт рекурсии мы можем воспользоваться операцией, которую мы определяем, на “хвостовом” ряде.
Результат деления $Q$ удовлетворяет соотношению:
$$F = Q * G$$
Переписав $F$, $G$ и $Q$ в нашем представлении, получим
| $f + x F_1$ | $=$ | $(q + x Q_1) * G = qG + x Q_1 * G = q(g + x G_1) + x Q_1 * G$ |
| $=$ | $q g + x (q G_1 + Q_1 * G)$ |
Следовательно
| $q$ | $=$ | $f/g$ |
| $Q_1$ | $=$ | $(F_1 - q G_1)/G$ |
Если $g = 0$ деление имеет смысл только в том случае, если и $f = 0$. Переведём на Haskell:
instance (Eq a, Fractional a) => Fractional (Ps a) where
(0 :+: fs) / (0 :+: gs) = fs / gs
(f :+: fs) / (g :+: gs) = q :+: ((fs - q .* gs)/(g :+: gs))
where q = f/g
fromRational x = p0 (fromRational x)
Производная одного члена ряда вычисляется так:
$$\frac{d}{dx} x^n = n x^{n-1}$$
Из этого выражения по свойствам производной
$$\frac{d}{dx}(f(x) + g(x)) = \frac{d}{dx}f(x) + \frac{d}{dx}g(x)$$
$$\frac{d}{dx} (k * f(x)) = k * \frac{d}{dx} f(x)$$
мы можем получить формулу для всего ряда:
$$\frac{d}{dx} F(x) = f_1 + 2 f_2 x + 3 f_3 x^2 + 4 f_4 x^3 + ...$$
Для реализации нам понадобится вспомогательная функция, которая будет обновлять значение дополнительного множителя $n$ в выражении $n x^{n-1}$:
diff :: Num a => Ps a -> Ps a
diff (f :+: fs) = diff' 1 fs
where diff' n (g :+: gs) = (n * g) :+: (diff' (n+1) gs)
Также мы можем вычислить и интеграл степенного ряда:
int :: Fractional a => Ps a -> Ps a
int (f :+: fs) = 0 :+: (int' 1 fs)
where int' n (g :+: gs) = (g / n) :+: (int' (n+1) gs)
Мы можем выразить элементарные функции через операции взятия производной и интегрирования. К примеру уравнение для $e^x$ выглядит так:
$$\frac{dy}{dx} = y$$
Проинтегрируем с начальным условием $y(0) = 1$:
$$y(x) = 1 + \int_0^x y(t) dt$$
Теперь переведём на Haskell:
expx = 1 + int expx
Кажется невероятным, но это и есть определение экспоненты. Так же мы можем определить и функции для синуса и косинуса:
| $\frac{d}{dx} \sin{x}$ | $=$ | $\cos{x},$ | $\sin(0) = 0,$ |
| $\frac{d}{dx} \cos{x}$ | $=$ | $- \sin{x},$ | $\cos(0) = 1$ |
Что приводит нас к:
sinx = int cosx cosx = 1 - int sinx
И это работает! Вычисление этих функций возможно за счёт того, что вне зависимости от аргумента функция int вернёт ряд, у которого первый коэффициент равен нулю. Это значение подхватывается и используется на следующем шаге рекурсивных вычислений.
Через синус и косинус мы можем определить тангенс:
tanx = sinx / cosx
В этом примере мы рассмотрим одну интересную технику рекурсивных вычислений, которая называется мемоизацией (memoization). Она заключается в том, что мы запоминаем все значения, с которыми вызывалась функция и, если с данным значением функция уже вычислялась, просто используем значение из памяти, а если значение ещё не вычислялось, вычисляем его и сохраняем.
В ленивых языках программирования для мемоизации функций часто используется такой приём. Мы сохраняем все значения функции в некотором контейнере, а затем обращаемся к элементам. При этом значения сохраняются в контейнере и не перевычисляются. Это происходит за счёт ленивых вычислений. Что интересно вычисляются не все значения, а лишь те, которые нам действительно нужны, те которые мы извлекаем из контейнера хотя бы один раз.
Посмотрим на такой классический пример. Вычисление чисел Фибоначчи. Каждое последующее число ряда Фибоначчи равно сумме двух предыдущих. Наивное определение выглядит так:
fib :: Int -> Int fib 0 = 0 fib 1 = 1 fib n = fib (n-1) + fib (n-2)
В этом определении число вычислений растёт экспоненциально. Для того чтобы вычислить fib n нам нужно вычислить fib (n-1) и fib (n-2), для того чтобы вычислить каждое из них нам нужно вычислить ещё два числа, и так вычисления удваиваются на каждом шаге. Если мы вызовем в интерпретаторе fib 40, то вычислитель зависнет. Что интересно в этой функции вычисления пересекаются, они могут быть переиспользованы. Например для вычисления fib (n-1) и fib (n-2) нужно вычислить fib (n-2) (снова), fib (n-3), fib (n-3) (снова) и fib (n-4).
Если мы сохраним все значения функции в списке, каждый вызов функции будет вычислен лишь один раз:
fib' :: Int -> Int
fib' n = fibs !! n
where fibs = 0 : 1 : zipWith (+) fibs (tail fibs)
Попробуем вычислить для 40:
*Fib> fib' 40 102334155 *Fib> fib' 4040 700852629
Вычисления происходят мгновенно. Если задача состоит из множества подзадач, которые самоподобны и для вычисления последующих подзадач используются решения из предыдущих, стоит задуматься об использовании мемоизации. Такие задачи называются задачами динамического программирования. Вычисление чисел Фибоначчи яркий пример задачи динамического программирования.
Рассмотрим такую задачу. Дана прямоугольная “карта местности”, в каждой клетке целым числом указана высота точки. Необходимо разметить местность по следующим правилам:
Из каждой клетки карты вода стекает не более чем в одном из четырёх возможных направлений (“север”, “юг”, “запад”, “восток”).
Если у клетки нет ни одного соседа с высотой меньше её собственной высоты, то эта клетка – водосток, и вода из неё никуда дальше не течёт.
Иначе вода из текущей клетки стекает на соседнюю клетку с минимальной высотой.
Если таких соседей несколько, то вода стекает по первому из возможных направлений из списка “на север”, “на запад”, “на восток”, “на юг”.
Все клетки из которых вода стекает в один и тот же водосток принадлежат к одному бассейну водосбора. Необходимо отметить на карте все бассейны. Решение этой задачи встретилось мне в статье Дмитрия Астапова “Рекурсия+мемоизация = динамическое программирование”. Здесь оно и приводится с незначительными изменениями.
Карта местности представлена в виде двумерного массива, в каждой клетке которого отмечена высота точки, нам необходимо получить двумерный массив того же размера, который вместо высот содержит метки водостоков. Мы будем отмечать их буквами латинского алфавита в том порядке, в котором они встречаются при обходе карты сверху вниз, слева направо. Например:
1 2 3 4 5 6 a a a b b b 7 8 9 2 4 5 a a b b b b 3 5 3 3 6 7 -> c c d b b e 6 4 5 5 3 1 f g d b e e 2 2 4 5 3 7 f g g h h e
Для представления двумерного массива мы воспользуемся типом Array из стандартного модуля Data.Array. Тип Array имеет два параметра:
data Array i a
Первый указывает на индекс, а второй на содержание. Массивы уже встречались нам в главе о типе ST. Напомню, что подразумевается, что этот тип является экземпляром класса Ix, который описывает целочисленные индексы, вспомним его определение:
class Ord a => Ix a where
range :: (a, a) -> [a]
index :: (a, a) -> a -> Int
inRange :: (a, a) -> a -> Bool
rangeSize :: (a, a) -> Int
Первый аргумент у всех этих функций это пара, которая представляет верхнюю и нижнюю грань последовательности. Попробуйте догадаться, что делают методы этого класса по типам и именам.
Для двумерного массива индекс будет задаваться парой целых чисел:
import Data.Array type Coord = (Int, Int) type HeightMap = Array Coord Int type SinkMap = Array Coord Coord type LabelMap = Array Coord Char
Значение типа HeightMap хранит карту высот, значение типа SinkMap хранит в каждой координате, ту точку, которая является водостоком для данной точки. Тип LabelMap обозначает итоговую разметку водостоков. Начнём с функции:
flow :: HeightMap -> SinkMap
Мы будем решать эту задачу рекурсивно. Представим, что мы знаем водостоки для всех точек кроме данной. Для каждой точки мы можем узнать в какую сторону из неё стекает вода. При этом водосток для следующей точки такой же как и для текущей. Если же из данной точки вода никуда не течёт, то она сама является водостоком. Мы определим эту функцию через комбинатор неподвижной точки fix.:
flow :: HeightMap -> SinkMap
flow arr = fix $ \result -> listArray (bounds arr) $
map (\x -> maybe x (result !) $ getSink arr x) $
range $ bounds arr
getSink :: HeightMap -> Coord -> Maybe Coord
Мы ищем решение в виде неподвижной точки функции, которая принимает карту стоков и возвращает карту стоков. Функция getSink по данной точке на карте вычисляет соседнюю точку, в которую стекает вода. Эта функция частично определена, поскольку для водостоков нет такой соседней точки, в которую бы утекала вода. Функция listArray конструирует значение типа Array из списка значений. Первым аргументом она принимает диапазон значений для индексов. Размеры массива совпадают с размерами карты высот, поэтому первым аргументом мы передаём bounds arr.
Теперь разберёмся с тем как заполняются значения в список. Сначала мы создаём список координат исходной карты высот с помощью выражения:
range $ bounds arr
После этого мы по координатам точек находим водостоки, причём сразу для всех точек. Это происходит в лямбда-функции:
\x -> maybe x (result !) $ getSink arr x
Мы принимаем текущую координату и с помощью функции getSink находим соседнюю точку, в которую убегает вода. Если такой точки нет, то в следующем выражении мы вернём исходную точку, поскольку в этом случае она и будет водостоком, а если такая соседняя точка всё-таки есть мы спросим результат из будущего. Мы обратимся к результату (result !), посмотрим каким окажется водосток для соседней точки и вернём это значение. Поскольку за счёт ленивых вычислений значения результирующего массива вычисляются лишь один раз, после того как мы найдём водосток для данной точки этим результатом смогут воспользоваться все соседние точки. При этом порядок обращения к значениям из будущих вычислений не играет роли.
Осталось только определить функцию поиска ближайшего стока и функцию разметки.
getSink :: HeightMap -> Coord -> Maybe Coord
getSink arr (x, y)
| null sinks = Nothing
| otherwise = Just $ snd $ minimum $ map (\i -> (arr!i, i)) sinks
where sinks = filter p [(x+1, y), (x-1, y), (x, y-1), (x, y+1)]
p i = inRange (bounds arr) i && arr ! i < arr ! (x, y)
В функции разметки мы воспользуемся ассоциативным массивом из модуля Data.Map. Функция nub из модуля Data.List убирает из списка повторяющиеся элементы. Затем мы составляем список пар из координат водостоков и меток и в самом конце размечаем исходный массив:
label :: SinkMap -> LabelMap
label a = fmap (m M.! ) a
where m = M.fromList $ flip zip ['a' .. ] $ nub $ elems a
Мы выяснили, что значение может редуцироваться только при сопоставлении с образцом и в специальной функции seq. Функцию seq мы можем применять, а можем и не применять. Но кажется, что в декомпозиции мы не можем уйти от необходимости проведения хотя бы одной редукции. Оказывается можем, в Haskell для этого предусмотрены специальные ленивые образцы (lazy patterns). Они обозначаются знаком тильда:
lazyHead :: [a] -> a lazyHead ~(x:xs) = x
Перед скобками сопоставления с образцом пишется символ тильда. Этим мы говорим вычислителю: доверься мне, здесь точно такой образец, можешь даже не проверять дальше. Он и правда дальше не пойдёт. Например если мы напишем такое определение:
lazySafeHead :: [a] -> Maybe a lazySafeHead ~(x:xs) = Just x lazySafeHead [] = Nothing
Если мы подставим в эту функцию пустой список мы получим ошибку времени выполнения, вычислитель доверился нам в первом уравнении, а мы его обманули. Сохраним в модуле Strict и проверим:
Prelude Strict> :! ghc --make Strict
[1 of 1] Compiling Strict ( Strict.hs, Strict.o )
Strict.hs:67:0:
Warning: Pattern match(es) are overlapped
In the definition of `lazySafeHead': lazySafeHead [] = ...
Prelude Strict> :l Strict
Ok, modules loaded: Strict.
Prelude Strict> lazySafeHead [1,2,3]
Just 1
Prelude Strict> lazySafeHead []
Just *** Exception: Strict.hs:(67,0)-(68,29): Irrefutable
pattern failed for pattern (x : xs)
При компиляции нам даже сообщили о том, что образцы в декомпозиции пересекаются. Но мы были упрямы и напоролись на ошибку, если мы поменяем образцы местами, то всё пройдёт гладко:
Prelude Strict> :! ghc --make Strict [1 of 1] Compiling Strict ( Strict.hs, Strict.o ) Prelude Strict> :l Strict Ok, modules loaded: Strict. Prelude Strict> lazySafeHead [] Nothing
Отметим, что сопоставление с образцом в let и where выражениях является ленивым. Функцию lazyHead мы могли бы написать и так:
lazyHead a = x
where (x:xs) = a
lazyHead a =
let (x:xs) = a
in x
Посмотрим как используются ленивые образцы при построении потоков, или бесконечных списков. Мы будем представлять функции одного аргумента потоками значений с одинаковым шагом. Так мы будем представлять непрерывные функции дискретными сигналами. Считаем, что шаг дискретизации (или шаг между соседними точками) нам известен.
$$f : R \rightarrow R \ \Rightarrow \ f_n = f(n \tau),\quad n = 0,1,2, ...$$
Где $\tau$ – шаг дискретизации, а $n$ пробегает все натуральные числа. Определим функцию решения дифференциальных уравнений вида:
$$\frac{dx}{dt} = f(t)$$
$$x(0) = \hat{x}$$
Символ $\hat{x}$ означает начальное значение функции $x$. Перейдём к дискретным сигналам:
$\frac{x_n - x_{n-1}}{\tau} = f_n, \quad x_0 = \hat{x}$
Где $\tau$ – шаг дискретизации, а $x$ и $f$ – это потоки чисел, индекс n пробегает от нуля до бесконечности по всем точкам функции, превращённой в дискретный сигнал. Такой метод приближения дифференциальных уравнений называют методом Эйлера. Теперь мы можем выразить следующий элемент сигнала через предыдущий.
$$x_n = x_{n-1} + \tau f_n, \quad x_0 = \hat{x}$$
Закодируем это уравнение:
-- шаг дискретизации dt :: Fractional a => a dt = 1e-3 -- метод Эйлера int :: Fractional a => a -> [a] -> [a] int x0 (f:fs) = x0 : int (x0 + dt * f) fs
Смотрите в функции int мы принимаем начальное значение x0 и поток всех значений функции правой части уравнения, поток значений функции $f(t)$. Мы помещаем начальное значение в первый элемент результата, а остальные значения получаем рекурсивно.
Определим две вспомогательные функции:
time :: (Enum a, Fractional a) => [a]
time = [0, dt ..]
dist :: Fractional a => Int -> [a] -> [a] -> a
dist n a b = ( / fromIntegral n) $
foldl' (+) 0 $ take n $ map abs $ zipWith (-) a b
Функция time пробегает все значения отсчётов шага дискретизации по времени. Это тождественная функция представленная в виде потока с шагом dt.
Функция проверки результата dist принимает два потока и по ним считает расстояние между ними. Эта функция говорит, что расстояние между двумя потоками в n первых точках равно сумме модулей разности между значениями потоков. Для того чтобы оценить среднее расхождение, мы делим в конце результат на число точек.
Также импортируем для удобства символьный синоним для fmap из модуля Control.Applicative.
import Control.Applicative((<$>)) ...
Проверим функцию int. Для этого сохраним все новые функции в модуле Stream.hs. Загрузим модуль в интерпретатор и вычислим производную какой-нибудь функции. Найдём решение для правой части константы и проверим, что у нас получилась тождественная функция:
*Stream> dist 1000 time $ int 0 $ repeat 1 7.37188088351104e-17
Функции практически совпадают, порядок ошибки составляет $10^{-16}$. Так и должно быть для линейных функций. Посмотрим, что будет если в правой части уравнения стоит тождественная функция:
*Stream> dist 1000 ((\t -> t^2/2) <$> time) $ int 0 time 2.497500000001403e-4
Решение этого уравнения равно функции $\frac{t^2}{2}$. Здесь мы видим, что результаты уже не такие хорошие.
Есть функции, которые определяются рекурсивно в терминах дифференциальных уравнений, например экспонента будет решением такого уравнения:
$$\frac{dx}{dt} = x$$
$$x(t) = x(0) + \int_{0}^{t} x(\tau) d \tau$$
Опишем это уравнение в Haskell:
e = int 1 e
Наше описание копирует исходное математическое определение. Добавим это уравнение в модуль Stream и проверим результаты:
*Stream> dist 1000 (map exp time) e ^CInterrupted.
К сожалению вычисление зависло. Нажмём ctrl+c и разберёмся почему. Для этого распишем вычисление потока чисел e:
e -- раскроем e
=> int 1 e -- раскроем int, во втором в аргументе
-- int стоит декомпозиция,
=> int 1 e@(f:fs) -- для того чтобы узнать какое уравнение
-- для int выбрать нам нужно раскрыть
-- второй аргумент, узнать корневой
-- конструктор, раскроем второй аргумент:
=> int 1 (int 1 e)
=> int 1 (int 1e@(f:fs)) -- такая же ситуация
=> int 1 (int 1 (int 1 e))
Проблема в том, что первый элемент решения мы знаем, мы передаём его первым аргументом и присоединяем к решению, но справа от знака равно. Но для того чтобы перейти в правую часть вычислителю нужно проверить все аргументы, в которых есть декомпозиция. И он начинает проверять, но слишком рано. Нам бы хотелось, чтобы он сначала присоединил к решению первый аргумент, а затем выполнял бы вычисления следующего элемента.
C помощью ленивых образцов мы можем отложить декомпозицию второго аргумента на потом:
int :: Fractional a => a -> [a] -> [a] int x0 ~(f:fs) = x0 : int (x0 + dt * f) fs
Теперь мы видим:
*Stream> dist 1000 (map exp time) e 4.988984990735441e-4
Вычисления происходят. С помощью взаимно-рекурсивных функций мы можем определить функции синус и косинус:
sinx = int 0 cosx cosx = int 1 (negate <$> sinx)
Эти функции описывают точку, которая бегает по окружности. Вот математическое определение:
| $\frac{dx}{dt}$ | $=$ | $y$ |
| $\frac{dy}{dt}$ | $=$ | $- x$ |
| $x(0)$ | $=$ | $0$ |
| $y(0)$ | $=$ | $1$ |
Проверим в интерпретаторе:
*Stream> dist 1000 (sin <$> time) sinx 1.5027460329809257e-4 *Stream> dist 1000 (cos <$> time) cosx 1.9088156807382827e-4
Так с помощью ленивых образцов нам удалось попасть в правую часть уравнения для функции int, не раскрывая до конца аргументы в левой части. С помощью этого мы могли ссылаться в сопоставлении с образцом на значение, которое ещё не было вычислено.
Ленивые вычисления повышают модульность программ. Мы можем в одной части программы создать все возможные решения, а в другой выбрать лучшие по какому-либо признаку. Также мы посмотрели на интересную технику написания рекурсивных функций, которая называется мемоизацией. Мемоизация означает, что мы не вычисляем повторно значения некоторой функции, а сохраняем их и используем в дальнейших вычислениях. Мы узнали новую синтаксическую конструкцию. Оказывается мы можем не только бороться с ленью, но и поощрять её. Лень поощряется ленивыми образцами. Они отменяют приведение к слабой заголовочной нормальной форме при декомпозиции аргументов. Они пишутся как обычные образцы, но со знаком тильда:
lazyHead ~(x:xs) = x
Мы говорим вычислителю: поверь мне, это значение может иметь только такой вид, потом посмотришь так ли это, когда значения тебе понадобятся. Поэтому ленивые образцы проходят сопоставление с образцом в любом случае.
Сопоставление с образцом в let и where выражениях является ленивым. Функцию lazyHead мы могли бы написать и так:
lazyHead a = x
where (x:xs) = a
lazyHead a =
let (x:xs) = a
in x
Мы побывали на выставке ленивых программ. Присмотритесь ещё раз к решениям задач этой главы и подумайте какую роль сыграли ленивые вычисления в каждом из случаев, какие мотивы обыгрываются в этих примерах. Также подумайте каким было бы решение, если бы в Haskell использовалась стратегия вычисления по значению. Критически настроенные читатели могут с помощью профилирования проверить эффективность программ из этой главы.