arma-thesis

git clone https://git.igankevich.com/arma-thesis.git
Log | Files | Refs | LICENSE

commit 18861e39fb2a100c40881430d94794a7863b3250
parent 99f19cefd2d4c3feb0f8525008daf6b650f579b5
Author: Ivan Gankevich <igankevich@ya.ru>
Date:   Wed, 15 Feb 2017 17:27:22 +0300

TODO Insert unbreakable spaces before ~cite~ and ~ref~, and ---.

Diffstat:
phd-diss-ru.org | 294+++++++++++++++++++++++++++++++++++++------------------------------------------
phd-diss.org | 212++++++++++++++++++++++++++++++++++++-------------------------------------------
2 files changed, 232 insertions(+), 274 deletions(-)

diff --git a/phd-diss-ru.org b/phd-diss-ru.org @@ -270,8 +270,7 @@ NaN: 29, -nan, 1.798e+36, -1.04284e+38, inf, -1.798e+36, -1.798e+36 для расчета качки судна, оценки величины воздействия внешних сил на плавучую платформу или другой морской объект, а также для оценки вероятности опрокидывания судна при заданных погодных условиях; однако, большинство из них -используют линейную теорию для моделирования морского волнения -cite:shin2003nonlinear,van2007forensic,kat2001prediction,van2002development, в +используют линейную теорию для моделирования морского волнения\nbsp{}cite:shin2003nonlinear,van2007forensic,kat2001prediction,van2002development, в рамках которой сложно воспроизвести определенные особенности ветроволнового климата. Среди них можно выделить переход от нормальных погодных условий к шторму и волнение, вызванное наложением множества систем ветровых волн и волн @@ -389,8 +388,7 @@ cite:shin2003nonlinear,van2007forensic,kat2001prediction,van2002development, в **** Методология и методы исследования. Программная реализация модели АРСС и формула вычисления давлений создавалась -поэтапно: прототип, написанный высокойровневом инженерном языке -cite:mathematica10,octave2015, был преобразован в программу на языке более +поэтапно: прототип, написанный высокойровневом инженерном языке\nbsp{}cite:mathematica10,octave2015, был преобразован в программу на языке более низкого уровня (C++). Реализация одних и тех же формул и алгоритмов на языках разного уровня (ввиду использования различных абстракций и языковых примитивов) позволяет выявить и исправить ошибки, которые остались бы незамеченными в случае @@ -436,20 +434,18 @@ Motion Programme (LAMP), программе для моделирования к **** Формула для поля давлений. Задача определения поля давлений под взволнованной морской поверхностью представляет собой обратную задачу гидродинамики для несжимаемой невязкой -жидкости. Система уравнений для нее в общем виде записывается как -cite:kochin1966theoretical +жидкости. Система уравнений для нее в общем виде записывается как\nbsp{}cite:kochin1966theoretical \begin{align} & \nabla^2\phi = 0,\nonumber\\ & \phi_t+\frac{1}{2} |\vec{\upsilon}|^2 + g\zeta=-\frac{p}{\rho}, & \text{на }z=\zeta(x,y,t),\label{eq:problem}\\ & D\zeta = \nabla \phi \cdot \vec{n}, & \text{на }z=\zeta(x,y,t),\nonumber \end{align} -где \(\phi\) --- потенциал скорости, \(\zeta\) --- подъем (аппликата) -взволнованной поверхности, \(p\) --- давление жидкости, \(\rho\) --- плотность -жидкости, \(\vec{\upsilon}=(\phi_x,\phi_y,\phi_z)\) --- вектор скорости, \(g\) ---- ускорение свободного падения и \(D\) --- субстанциональная производная +где \(\phi\)\nbsp{}--- потенциал скорости, \(\zeta\)\nbsp{}--- подъем (аппликата) +взволнованной поверхности, \(p\)\nbsp{}--- давление жидкости, \(\rho\)\nbsp{}--- плотность +жидкости, \(\vec{\upsilon}=(\phi_x,\phi_y,\phi_z)\)\nbsp{}--- вектор скорости, \(g\)\nbsp{}--- ускорение свободного падения и \(D\)\nbsp{}--- субстанциональная производная (производная Лагранжа). Первое уравнение является уравнением неразрывности -(уравнение Лапласа), второе --- законом сохранения импульса (которое иногда -называют динамическим граничным условием); третье уравнение --- кинематическое +(уравнение Лапласа), второе\nbsp{}--- законом сохранения импульса (которое иногда +называют динамическим граничным условием); третье уравнение\nbsp{}--- кинематическое граничное условие, которое сводится к равенству скорости перемещения этой поверхности (\(D\zeta\)) нормальной составляющей скорости жидкости (\(\nabla\phi\cdot\vec{n}\)). @@ -459,7 +455,7 @@ cite:kochin1966theoretical формулой для определения поля давлений по значениям производных потенциалов скорости, полученных из оставшихся уравнений. Таким образом, с математической точки зрения обратная задача гидродинамики сводится к решению уравнения Лапласа -со смешанным ГУ --- задаче Робена. +со смешанным ГУ\nbsp{}--- задаче Робена. * Обзор литературы ** Анализ моделей морского волнения @@ -472,9 +468,8 @@ cite:kochin1966theoretical *** Модель Лонге---Хиггинса Наиболее простой моделью, формула которой выводится в рамках линейной теории -волн, является модель Лонге---Хиггинса (ЛХ) cite:longuet1957statistical. -Подробный сравнительный анализ этой модели и модели АРСС проведен в работах -cite:degtyarev2011modelling,boukhanovsky1997thesis. +волн, является модель Лонге---Хиггинса (ЛХ)\nbsp{}cite:longuet1957statistical. +Подробный сравнительный анализ этой модели и модели АРСС проведен в работах\nbsp{}cite:degtyarev2011modelling,boukhanovsky1997thesis. Модель ЛХ представляет взволнованную морскую поверхность в виде суперпозиции элементарных гармонических волн случайных амплитуд \(c_n\) и фаз \(\epsilon_n\), @@ -493,7 +488,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. \begin{equation*} 2E_\zeta(u,v)\, du\, dv = \sum\limits_n c_n^2, \end{equation*} -где \(E_\zeta(u,v)\) --- двумерная спектральная плотность энергии волн. +где \(E_\zeta(u,v)\)\nbsp{}--- двумерная спектральная плотность энергии волн. Коэффициенты \(c_n\) определяются из энергетического спектра волнения \(S(\omega)\) по формуле \begin{equation*} @@ -510,8 +505,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. независимости от спектра, подаваемого на вход модели. Использование меньшего количества коэффициентов может решить проблему, но также уменьшит период реализации. Таким образом, использование модели ЛХ для генерации волн с - негауссовым распределением аппликат (которое имеют реальные морские волны - cite:huang1980experimental,рожков1996теория) не реализуемо на практике. + негауссовым распределением аппликат (которое имеют реальные морские волны\nbsp{}cite:huang1980experimental,рожков1996теория) не реализуемо на практике. 2. С вычислительной точки зрения, недостатком модели является нелинейный рост времени генерации поверхности с увеличением размера реализации. Чем больше размер реализации, тем больше коэффициентов (дискретных точек @@ -524,7 +518,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. - В программной реализации скорость сходимости выражения ([[eq:longuet-higgins]]) может быть низкой, т.к. фазы \(\epsilon_n\) имеют вероятностный характер. - Обобщение модели для негауссовых и нелинейных процессов сопряжено с большой - трудоемкостью вычислений cite:рожков1990вероятностные. + трудоемкостью вычислений\nbsp{}cite:рожков1990вероятностные. Таким образом, модель Лонге---Хиггинса применима для решения задачи генерации взволнованной морской поверхности только в линейной постановке (в рамках теории @@ -533,7 +527,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. более совершенных моделей. *** Модель АРСС -В cite:spanos1982arma модель АРСС используется для генерации временного ряда, +В\nbsp{}cite:spanos1982arma модель АРСС используется для генерации временного ряда, спектр которого совпадает с аппроксимацией Пирсона---Московица для спектров морского волнения. Авторы проводят эксперименты для одномерных моделей АР, СС и АРСС. Они отмечают превосходное совпадение полученного и исходного спектров и @@ -541,7 +535,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. моделями, основанными на суммировании большого числа гармоник со случайными фазами. Также отмечается, что для того чтобы спектр полученного временного ряда совпадал с заданным, модели СС требуется меньшее количество коэффициентов, чем -модели АР. В cite:spanos1996efficient автор обобщает формулы для нахождения +модели АР. В\nbsp{}cite:spanos1996efficient автор обобщает формулы для нахождения коэффициентов модели АРСС для случая нескольких (векторов) переменных. Отличие данной работы от вышеперечисленных отличается в исследовании трехмерной @@ -565,7 +559,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. а не является абстрактным многомерным случайным процессом, совпадающим с реальным лишь статистически. -В cite:fusco2010short модель АР используется для прогнозирования волн зыби для +В\nbsp{}cite:fusco2010short модель АР используется для прогнозирования волн зыби для управления преобразователем энергии волн (ПЭВ) в реальном времени. Для эффективной работы ПЭВ необходимо чтобы частота встроенного осциллятора совпадала с частотой морских волн. Авторы статьи представляют подъем волны как @@ -576,7 +570,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. ** Известные формулы определения поля давлений *** Теория волн малых амплитуд -В cite:stab2012,детярев1998моделирование,degtyarev1997analysis дается решение +В\nbsp{}cite:stab2012,детярев1998моделирование,degtyarev1997analysis дается решение обратной задачи гидродинамики для случая идеальной несжимаемой жидкости в рамках теории волн малых амплитуд (в предположении, что длина волны много больше ее высоты: \(\lambda \gg h\)). В этом случае обратная задача линейна и сводится к @@ -585,7 +579,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. потенциала скорости. Предположение о малости амплитуд волн означает слабое изменение локального волнового числа во времени и пространстве по сравнению с подъемом (аппликатой) взволнованной поверхности. Это позволяет вычислить -производную подъема поверхности по \(z\) как \(\zeta_z=k\zeta\), где \(k\) --- +производную подъема поверхности по \(z\) как \(\zeta_z=k\zeta\), где \(k\)\nbsp{}--- волновое число. В двухмерном случае решение записывается явной формулой \begin{align} \left.\frac{\partial\phi}{\partial x}\right|_{x,t}= & @@ -594,7 +588,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. z+\alpha\dot{\alpha}}{\sqrt{1+\alpha^{2}}}e^{I(x)}dx,\label{eq:old-sol-2d}\\ I(x)= & \int\limits_{0}^x\frac{\partial\alpha/\partial z}{1+\alpha^{2}}dx,\nonumber \end{align} -где \(\alpha\) --- уклоны волн. В трехмерном случае решение записывается в виде +где \(\alpha\)\nbsp{}--- уклоны волн. В трехмерном случае решение записывается в виде эллиптического дифференциального уравнения в частных производных \begin{align*} & \frac{\partial^2 \phi}{\partial x^2} \left( 1 + \alpha_x^2 \right) + @@ -656,10 +650,10 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. , \label{eq:arma-process} \end{equation} -где \(\zeta\) --- подъем (аппликата) взволнованной поверхности, \(\Phi\) --- -коэффициенты процесса АР, \(\Theta\) --- коэффициенты процесса СС, \(\epsilon\) --- -белый шум, имеющий Гауссово распределение, \(\vec N\) --- порядок процесса АР, -\(\vec M\) --- порядок процесса СС, причем \(\Phi_{\vec{0}}\equiv0\), +где \(\zeta\)\nbsp{}--- подъем (аппликата) взволнованной поверхности, \(\Phi\)\nbsp{}--- +коэффициенты процесса АР, \(\Theta\)\nbsp{}--- коэффициенты процесса СС, \(\epsilon\)\nbsp{}--- +белый шум, имеющий Гауссово распределение, \(\vec N\)\nbsp{}--- порядок процесса АР, +\(\vec M\)\nbsp{}--- порядок процесса СС, причем \(\Phi_{\vec{0}}\equiv0\), \(\Theta_{\vec{0}}\equiv0\). Здесь стрелки обозначают многокомпонентные индексы, содержащие отдельную компоненту для каждого измерения. В общем случае в качестве компонент могут выступать любые скалярные величины (температура, соленость, @@ -667,7 +661,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. коэффициенты и порядки процессов АР и СС. **** Процесс авторегрессии (АР). -Процесс АР --- это процесс АРСС только лишь с одним случайным импульсом вместо их +Процесс АР\nbsp{}--- это процесс АРСС только лишь с одним случайным импульсом вместо их взвешенной суммы: \begin{equation} \zeta_{\vec i} @@ -699,7 +693,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. 0, \quad \text{if } \vec{k}\neq0, \end{cases} \end{equation} -где \(\gamma\) --- АКФ процесса \(\zeta\), \(\Var{\epsilon}\) --- дисперсия +где \(\gamma\)\nbsp{}--- АКФ процесса \(\zeta\), \(\Var{\epsilon}\)\nbsp{}--- дисперсия белого шума. Матричная форма трехмерной системы уравнений Юла---Уокера, используемой в данной работе, имеет следующий вид. \begin{equation*} @@ -772,7 +766,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. \end{equation*} **** Процесс скользящего среднего (СС). -Процесс СС --- это процесс АРСС, в котором \(\Phi\equiv0\): +Процесс СС\nbsp{}--- это процесс АРСС, в котором \(\Phi\equiv0\): \begin{equation} \zeta_{\vec i} = @@ -809,7 +803,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. \Theta_{\vec j}^2 }. \end{equation*} -Авторы cite:box1976time предлагают использовать метод Ньютона---Рафсона для +Авторы\nbsp{}cite:box1976time предлагают использовать метод Ньютона---Рафсона для решения этого уравнения с большей точностью, однако, этот метод не подходит для трех измерений. Использование более медленного метода не оказывает большого эффекта на общую производительность программы, потому что количество @@ -826,7 +820,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. АР, однако, для того чтобы АКФ результирующего процесса соответствовала заданной, необходимо предварительно скорректировать значения коэффициентов АР. Существует несколько способов "смешивания" процессов АР и СС. -- Подход, предложенный авторами cite:box1976time, который включается в себя +- Подход, предложенный авторами\nbsp{}cite:box1976time, который включается в себя разделение АКФ на часть для процесса АР и часть для процесса СС по каждому из измерений, не подходит в данной ситуации, поскольку в трех измерениях невозможно таким образом разделить АКФ: всегда останутся части, которые не @@ -835,7 +829,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. АКФ для обоих процессов разных порядков, однако, тогда характеристики реализации (математической ожидание, дисперсия и др.) будут смещены: они станут характеристика двух наложенных друг на друга процессов. -Для первого подхода авторами cite:box1976time предложена формула корректировки +Для первого подхода авторами\nbsp{}cite:box1976time предложена формула корректировки коэффициентов процесса АР, для второго же подхода такой формулы нет. Таким образом, лучшим решением на данный момент является использование процессов АР и СС по отдельности. @@ -844,22 +838,22 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. Одной из проблем в применении модели АРСС для генерации взволнованной морской поверхности является то, что для разных профилей волн /необходимо/ использовать разные процессы: стоячие волны моделируются только процессом АР, а прогрессивные -волны --- только процессом СС. Это утверждение пришло из практики: если +волны\nbsp{}--- только процессом СС. Это утверждение пришло из практики: если попытаться использовать процессы наоборот, результирующая реализация либо расходится, либо не представляет собой реальные морские волны (такое происходит в случае необратимого процесса СС, который всегда стационарен). Таким образом, процесс АР может быть использован только для моделирования стоячих волн, а -процесс СС --- для прогрессивных волн. +процесс СС\nbsp{}--- для прогрессивных волн. Другой проблемой является невозможность автоматического определения оптимального количества коэффициентов для трехмерных процессов АР и СС. Для одномерных -процессов существуют итеративные методы cite:box1976time, однако они расходятся +процессов существуют итеративные методы\nbsp{}cite:box1976time, однако они расходятся в трехмерном случае. Последней проблемой, которая описана в разделе [[#sec:how-to-mix-ARMA]], является невозможность "смешать" процесс АР и СС в трех измерениях. -Практика показывает, что некоторые утверждения авторов cite:box1976time не +Практика показывает, что некоторые утверждения авторов\nbsp{}cite:box1976time не выполняются для трехмерной модели АРСС. Например, авторы утверждают, что АКФ процесса СС обрывается на отсчете \(q\), а АКФ процесса АР затухает на бесконечности, однако, на практике при использовании слабо затухающей и @@ -883,18 +877,18 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. Модель АРСС позволяет учесть асимметричность распределения волновых аппликат, т.е. генерировать морские волны, закон распределения аппликат которых имеет ненулевой эксцесс и асимметрию. Такой закон распределения характерен для реальных -морских волн cite:longuet1963nonlinear. +морских волн\nbsp{}cite:longuet1963nonlinear. Асимметричность волн моделируется с помощью нелинейного безынерционного преобразования (НБП) случайного процесса, однако, любое нелинейное преобразование случайного процесса приводит к преобразованию его АКФ. Для того чтобы подавить этот эффект, необходимо предварительно преобразовать АКФ, как -показано в cite:boukhanovsky1997thesis. +показано в\nbsp{}cite:boukhanovsky1997thesis. **** Преобразование взволнованной поверхности. Формула \(z=f(y)\) преобразования взволнованной поверхности к необходимому одномерному закону распределения \(F(z)\) получается путем решения нелинейного -трансцендентного уравнения \(F(z) = \Phi(y)\), где \(\Phi(y)\) --- функция +трансцендентного уравнения \(F(z) = \Phi(y)\), где \(\Phi(y)\)\nbsp{}--- функция одномерного нормального закона распределения. Поскольку функция распределения аппликат морских волн часто задается некоторой аппроксимацией, основанной на натурных данных, то это уравнение целесообразно решать численно в каждой точке @@ -926,7 +920,7 @@ cite:degtyarev2011modelling,boukhanovsky1997thesis. \int\limits_{0}^\infty f(y) H_m(y) \exp\left[ -\frac{y^2}{2} \right], \end{equation*} -\(H_m\) --- полином Эрмита, а \(f(y)\) --- решение уравнения +\(H_m\)\nbsp{}--- полином Эрмита, а \(f(y)\)\nbsp{}--- решение уравнения eqref:eq:distribution-transformation. Воспользовавшись полиномиальной аппроксимацией \(f(y) \approx \sum\limits_i d_i y^i\) и аналитическими выражениями для полнимов Эрмита, формулу определения коэффициентов можно упростить, @@ -949,7 +943,7 @@ eqref:eq:distribution-transformation. Воспользовавшись поли \frac{C_k^2}{k!} \right| \leq \epsilon. \end{equation*} -В cite:boukhanovsky1997thesis автор предлагает использовать полиномиальную +В\nbsp{}cite:boukhanovsky1997thesis автор предлагает использовать полиномиальную аппроксимацию для \(f(y)\) также для преобразования поверхности, однако на практике в реализации взволнованной поверхности часто находятся точки, выпадающие за промежуток на котором построена аппроксимация, что приводит к @@ -984,8 +978,8 @@ eqref:eq:distribution-transformation эффективнее решать мет частных производных преобразуется в неявную разностную схему, решаемую итерационным методом, на каждом шаге которого ищется решение трехдиагональной или пятидиагональной СЛАУ методом прогонки (алгоритм Томаса). Асимптотическая -сложность алгоритма составляет \(\mathcal{O}({n}{m})\), где \(n\) --- количество -точек на сетке взволнованной поверхности, \(m\) --- число итераций. Несмотря на +сложность алгоритма составляет \(\mathcal{O}({n}{m})\), где \(n\)\nbsp{}--- количество +точек на сетке взволнованной поверхности, \(m\)\nbsp{}--- число итераций. Несмотря на широкое распространение, итеративные алгоритмы неэффективно отображаются на архитектуру параллельных машин; в частности, отображение на сопроцессоры может включать в себя копирование данных на сопроцессор и обратно на каждой итерации, @@ -1045,7 +1039,7 @@ eqref:eq:distribution-transformation эффективнее решать мет \ast \InverseFourierY{E(u)}{x}, \end{equation*} -где \(\Fun{z}\) --- некоторая функция, вид которой будет определен в +где \(\Fun{z}\)\nbsp{}--- некоторая функция, вид которой будет определен в [[#sec:compute-delta]] и для которой выполняется соотношение \(\FourierY{\Fun{z}}{u}=e^{2\pi{u}{z}}\). Подставляя выражение для \(\phi\) в граничное условие, получим @@ -1100,7 +1094,7 @@ eqref:eq:guessed-sol-2d, получаем окончательное выраж **** Формула для жидкости конечной глубины. На дне водоема вертикальная составляющая скорости перемещения жидкости должна -равняться нулю, т.е. \(\phi_z=0\) на \(z=-h\), где \(h\) --- глубина водоема. В этом +равняться нулю, т.е. \(\phi_z=0\) на \(z=-h\), где \(h\)\nbsp{}--- глубина водоема. В этом случае пренебречь равенством \(v = -i u\), полученным из уравнения Лапласа, нельзя, и решение ищется в виде \begin{equation} @@ -1149,15 +1143,14 @@ eqref:eq:guessed-sol-2d, получаем окончательное выраж } \label{eq:solution-2d-full} \end{equation} -где \(\FunSecond{z}\) --- некоторая функция, вид которой будет определен в +где \(\FunSecond{z}\)\nbsp{}--- некоторая функция, вид которой будет определен в [[#sec:compute-delta]] и для которой выполняется соотношение \(\FourierY{\FunSecond{z}}{u}=\Sinh{2\pi{u}{z}}\). **** Сведение к формулам линейной теории волн. Справедливость полученных формул проверим, подставив в качестве \(\zeta(x,t)\) известные аналитические выражения для плоских волн. Символьные вычисления -преобразований Фурье в этом разделе производились с помощью пакета Mathematica -cite:mathematica10. В линейной теории широко используется предположение о +преобразований Фурье в этом разделе производились с помощью пакета Mathematica\nbsp{}cite:mathematica10. В линейной теории широко используется предположение о малости амплитуд волн, что позволяет упростить исходную систему уравнений eqref:eq:problem-2d до \begin{align*} @@ -1304,14 +1297,14 @@ eqref:eq:solution-2d-full до ее лучше всего вычислять с помощью программы для символьных вычислений. Затем, для практического применения она может быть аппроксимирована суперпозицией экспоненциально затухающих косинусов (именно так выглядит АКФ стационарного -процесса АРСС cite:box1976time). +процесса АРСС\nbsp{}cite:box1976time). **** Эмпирический метод. Впрочем, для трехмерного случая существует более простой эмпирический метод нахождения формы АКФ, не требующий использования сложного программного обеспечения. Известно, что АКФ, представляющая собой суперпозицию экспоненциально затухающих косинусов, является решением уравнения Стокса для -гравитационных волн cite:boccotti1983wind. Значит, если в моделируемом морском +гравитационных волн\nbsp{}cite:boccotti1983wind. Значит, если в моделируемом морском волнении важна только форма волны, а не точные ее характеристики, то заданный волновой профиль можно просто домножить на затухающую экспоненту, чтобы получить подходящую АКФ. Эта АКФ не отражает параметры волн, такие как высота и период, @@ -1408,9 +1401,9 @@ eqref:eq:solution-2d-full до полиномиальной аппроксимацией натурных данных, либо аналитически. **** Разложение в ряд Грама---Шарлье. -В cite:huang1980experimental было экспериментально показано, что распределение +В\nbsp{}cite:huang1980experimental было экспериментально показано, что распределение аппликат морской поверхности отличается от нормального ненулевым эксцессом и -асимметрией. В cite:рожков1996теория показано, что такое распределение +асимметрией. В\nbsp{}cite:рожков1996теория показано, что такое распределение раскладывается в ряд Грама---Шарлье: \begin{align} \label{eq:skew-normal-1} @@ -1434,12 +1427,12 @@ eqref:eq:solution-2d-full до +1 \right], \end{align} -где \(\phi(z)=\frac{1}{2}\mathrm{erf}(z/\sqrt{2})\), \(\gamma_1\) --- асимметрия, -\(\gamma_2\) --- эксцесс, \(f\) --- ФПР, \(F\) --- функция распределения (ФР). -Согласно cite:рожков1990вероятностные для аппликат морских волн значение +где \(\phi(z)=\frac{1}{2}\mathrm{erf}(z/\sqrt{2})\), \(\gamma_1\)\nbsp{}--- асимметрия, +\(\gamma_2\)\nbsp{}--- эксцесс, \(f\)\nbsp{}--- ФПР, \(F\)\nbsp{}--- функция распределения (ФР). +Согласно\nbsp{}cite:рожков1990вероятностные для аппликат морских волн значение асимметрии выбирается на интервале \(0,1\leq\gamma_1\leq{0,52}]\), а значение эксцесса на интервале \(0,1\leq\gamma_2\leq{0,7}\). Семейство плотностей -распределения при различных параметрах показано на [[fig:skew-normal-1]]. +распределения при различных параметрах показано на рис.\nbsp{}[[fig:skew-normal-1]]. #+name: fig:skew-normal-1 #+begin_src R :results output graphics :exports results :file build/skew-normal-1-ru.pdf @@ -1481,12 +1474,12 @@ legend( f(z; \alpha) & = \frac{e^{-\frac{z^2}{2}}}{\sqrt{2 \pi }} \mathrm{erfc}\left[-\frac{\alpha z}{\sqrt{2}}\right], \end{align} -где \(T\) --- функция Оуэна cite:owen1956tables. Эта формула не позволяет задать -значения асимметрии и эксцесса по отдельности --- оба значения регулируются +где \(T\)\nbsp{}--- функция Оуэна\nbsp{}cite:owen1956tables. Эта формула не позволяет задать +значения асимметрии и эксцесса по отдельности\nbsp{}--- оба значения регулируются параметром \(\alpha\). Преимущество данной формулы лишь в относительной простоте вычисления: эта функция встроена в некоторые программы и библиотеки математических функций. График функции для разных значений \(\alpha\) представлен -на [[fig:skew-normal-2]]. +на рис.\nbsp{}[[fig:skew-normal-2]]. #+name: fig:skew-normal-2 #+begin_src R :results output graphics :exports results :file build/skew-normal-2-ru.pdf @@ -1540,15 +1533,14 @@ legend( Чтобы исключить периодичность из сгенерированной моделью ветрового волнения реализации взволнованной поверхности, для генерации белого шума нужно использовать ГПСЧ с достаточно большим периодом. В качестве такого генератора в -работе используется параллельная реализация вихря Мерсенна -cite:matsumoto1998mersenne с периодом \(2^{19937}-1\). Это позволяет создавать +работе используется параллельная реализация вихря Мерсенна\nbsp{}cite:matsumoto1998mersenne с периодом \(2^{19937}-1\). Это позволяет создавать апериодичные реализации взволнованной морской поверхности для любых сценариев применения, встречаемых на практике. Запуск нескольких ГПСЧ с разными начальными состояниями в параллельных потоках не гарантирует некоррелированность генерируемых последовательностей псевдослучайных чисел, однако, можно воспользоваться алгоритмом динамического -создания вихрей Мерсенна cite:matsumoto1998dynamic, чтобы дать такую гарантию. +создания вихрей Мерсенна\nbsp{}cite:matsumoto1998dynamic, чтобы дать такую гарантию. Суть алгоритма заключается в поиске таких матриц начальных состояний генераторов, которые бы дали максимально некоррелированные последовательности псевдослучайных чисел при параллельном запуске нескольких вихрей Мерсенна с @@ -1561,7 +1553,7 @@ cite:matsumoto1998mersenne с периодом \(2^{19937}-1\). Это позв *** Алгоритм генерации взволнованной поверхности В модели АРСС значение подъема взволнованной поверхности в каждой точке зависит от предыдущих по пространству и времени значений, из-за чего в начале реализации -образуется так называемый /интервал разгона/ (см. рис. [[fig:ramp-up-interval]]) --- +образуется так называемый /интервал разгона/ (см. рис.\nbsp{}[[fig:ramp-up-interval]])\nbsp{}--- промежуток, на котором реализация не соответствует заданной АКФ. Способ решения этой проблемы зависит от контекста, в котором происходит моделирование. @@ -1577,10 +1569,9 @@ cite:matsumoto1998mersenne с периодом \(2^{19937}-1\). Это позв В алгоритме генерации взволнованной поверхности используется параллелизм по данным: реализация делится на равные части, каждая из которых генерируется -независимо, --- однако, в начале каждой из частей также присутствует интервал +независимо,\nbsp{}--- однако, в начале каждой из частей также присутствует интервал разгона. Для его исключения используется метод /сшивания/, часто применяемый в -обработке цифровых сигналов -cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms. Суть метода +обработке цифровых сигналов\nbsp{}cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms. Суть метода заключается в добавлении интервала равного по размеру интервалу разгона в конец каждой из частей. Затем взволнованная поверхность генерируется в каждой точки каждой из частей (включая добавленный интервал), интервал в конце части \(N\) @@ -1607,7 +1598,7 @@ arma.plot_ramp_up_interval(label="Интервал разгона") \(\Fun{z}=\InverseFourierY{e^{2\pi{u}{z}}}{x}\) и \(\FunSecond{z}=\InverseFourierY{\Sinh{2\pi{u}{z}}}{x}\), которые могут быть записаны аналитически различными выражениями и представляют сложность при -вычислении на компьютере. Каждая функция --- это преобразование Фурье от +вычислении на компьютере. Каждая функция\nbsp{}--- это преобразование Фурье от линейной комбинации экспонент, которое сводится к плохо определенной дельта функции комплексного аргумента (см. [[tab:delta-functions]]). Обычно такого типа функции записывают как произведение дельта функций от действительной и мнимой @@ -1633,8 +1624,7 @@ arma.plot_ramp_up_interval(label="Интервал разгона") :CUSTOM_ID: sec:verification :END: -Для модели АР в работах -cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis +Для модели АР в работах\nbsp{}cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis экспериментальным путем были верифицированы - распределения различных характеристик волн (высоты волн, длины волн, длины гребней, период волн, уклон волн, показатель трехмерности), @@ -1644,12 +1634,12 @@ cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis распределений различных характеристик волн. *** Верификация интегральных характеристик взволнованной поверхности -В cite:рожков1990вероятностные авторы показывают, что некоторые характеристики +В\nbsp{}cite:рожков1990вероятностные авторы показывают, что некоторые характеристики морских волн (перечисленные в таблице [[tab:weibull-shape]]) имеют распределение -Вейбулла, а подъем взволнованной поверхности --- нормальное распределение. Для +Вейбулла, а подъем взволнованной поверхности\nbsp{}--- нормальное распределение. Для верификации генерируемых моделями АР и СС реализаций используются спрямленные диаграммы (графики, в которых по оси \(OX\) откладываются квантили функции -распределения, вычисленные аналитически, а по оси \(OY\) --- вычисленные +распределения, вычисленные аналитически, а по оси \(OY\)\nbsp{}--- вычисленные экспериментально). Если экспериментально полученное распределение соответствует аналитическому, то график представляет собой прямую линию. Концы графика могут отклоняться от прямой линии, поскольку не могут быть надежно получены из @@ -1741,8 +1731,7 @@ eqref:eq:solution-2d-full с известными формулами линей **** Отличие от формул линейной теории волн. Эксперимент показывает, что поля потенциалов скоростей, полученные по формуле eqref:eq:solution-2d-full для конечной глубины и по формуле -eqref:eq:solution-2d-linear линейной теории, качественно отличаются (см. -[[fig:potential-field-nonlinear]]). Во-первых, контуры потенциала скорости имеют вид +eqref:eq:solution-2d-linear линейной теории, качественно отличаются (см. рис.\nbsp{}[[fig:potential-field-nonlinear]]). Во-первых, контуры потенциала скорости имеют вид затухающей синусоиды, что отличается от овальной формы, описываемой линейной теории волн. Во-вторых, по мере приближения к дну водоема потенциал скорости затухает гораздо быстрее, чем в линейной теории, а область, где сконцентрирована @@ -1765,13 +1754,12 @@ eqref:eq:solution-2d-linear линейной теории, качественн eqref:eq:solution-2d-full и формуле для волн малой амплитуды eqref:eq:old-sol-2d, сопоставимы для волн малых амплитуд. В этом эксперименте используются две реализации взволнованной морской поверхности, полученные по -модели АР: одна содержит волны малой амплитуды, другая --- большой. +модели АР: одна содержит волны малой амплитуды, другая\nbsp{}--- большой. Интегрирование в формуле eqref:eq:solution-2d-full ведется диапазону волновых чисел, полученному из морской поверхности. Для волн малой амплитуды обе формулы показывают сопоставимые результаты (разница в значениях скорости приписывается стохастической природе модели АР), в то время как для волн больших амплитуд -устойчивое поле скоростей дает только формула eqref:eq:solution-2d-full (рис. -[[fig:velocity-field-2d]]). Таким образом, общая формула eqref:eq:solution-2d-full +устойчивое поле скоростей дает только формула eqref:eq:solution-2d-full (рис. рис.\nbsp{}[[fig:velocity-field-2d]]). Таким образом, общая формула eqref:eq:solution-2d-full показывает удовлетворительные результаты, не вводя ограничения на амплитуду волн. @@ -1802,11 +1790,10 @@ eqref:eq:old-sol-2d, сопоставимы для волн малых ампл выходным данным предыдущего звена. Звенья конвейера распределяются по узлам вычислительного кластера, чтобы сделать возможным параллелизм по операциям, а затем данные, перемещающиеся между звеньями конвейера распределяются между -ядрами процессора, чтобы сделать возможным параллелизм по данным. На -[[fig:pipeline]] представлена схема конвейера обработки данных, в которой +ядрами процессора, чтобы сделать возможным параллелизм по данным. На рис.\nbsp{}[[fig:pipeline]] представлена схема конвейера обработки данных, в которой прямоугольниками со скругленными углами обозначены звенья конвейера, обычными -прямоугольниками --- массивы объектов из предметной области задачи, передаваемые -от одного звена к другому, а стрелками --- направление передачи данных. +прямоугольниками\nbsp{}--- массивы объектов из предметной области задачи, передаваемые +от одного звена к другому, а стрелками\nbsp{}--- направление передачи данных. Некоторые звенья разделены на /секции/, каждая из которых обрабатывает отдельную часть массива. Если звенья соединены без использования /барьера/ (горизонтальная или вертикальная полоса), то передача отдельных объектов между такими звеньями @@ -1925,9 +1912,7 @@ digraph { #+RESULTS: fig:pipeline [[file:build/pipeline-ru.pdf]] -Конвейер объектов можно считать развитием модели BSP (Bulk Synchronous Parallel) -cite:valiant1990bridging, применяемой в системах обработки графов -cite:malewicz2010pregel,seo2010hama. Конвейер позволяет исключить глобальную +Конвейер объектов можно считать развитием модели BSP (Bulk Synchronous Parallel)\nbsp{}cite:valiant1990bridging, применяемой в системах обработки графов\nbsp{}cite:malewicz2010pregel,seo2010hama. Конвейер позволяет исключить глобальную синхронизацию (где это возможно) между последовательно идущим этапами вычислений путем передачи данных между звеньев параллельно с вычислениями, в то время как в модели BSP глобальная синхронизация происходит после каждого шага. @@ -1982,8 +1967,7 @@ cite:malewicz2010pregel,seo2010hama. Конвейер позволяет иск должна быть ответственна за повторное выполнение не удавшейся задачи на одном из выживших узлов, тривиально. Чтобы повторно выполнить задачу на вершине иерархии, создается резервная задача, выполняющаяся на другом узле. Существует ряд систем, -которые способны выполнять направленные ациклические графы задач параллельно -cite:acun2014charmpp,islam2012oozie, но графы не подходят для определения +которые способны выполнять направленные ациклические графы задач параллельно\nbsp{}cite:acun2014charmpp,islam2012oozie, но графы не подходят для определения отношений руководитель-подчиненный между задачами, поскольку узел графа может иметь несколько родительских узлов. @@ -1993,10 +1977,10 @@ cite:acun2014charmpp,islam2012oozie, но графы не подходят дл к поломкам оборудования, т.е. обеспечение отказоустойчивости и высокой доступности, которое прозрачно для программиста. Реализация модели состоит из двух слоев: на нижнем слое находятся подпрограммы и классы для приложений, -работающих на одном узле (без сетевых взаимодействий), на верхнем слое --- для +работающих на одном узле (без сетевых взаимодействий), на верхнем слое\nbsp{}--- для приложений, работающих на произвольном количестве узлов. Модель включает в себя -два вида сильно связанных друг с другом сущностей --- /управляющие объекты/ (или -/ядра/) и /конвейеры/, --- которые используются совместно для написания +два вида сильно связанных друг с другом сущностей\nbsp{}--- /управляющие объекты/ (или +/ядра/) и /конвейеры/,\nbsp{}--- которые используются совместно для написания программы. Управляющие объекты реализуют логику (порядок выполнения) программы в методах @@ -2023,15 +2007,15 @@ cite:acun2014charmpp,islam2012oozie, но графы не подходят дл параграфе. Для каждого устройства используется отдельный конвейер. Существуют конвейеры для параллельной обработки, обработки по расписанию (периодические и отложенные задачи) и промежуточный конвейер для обработки управляющих объектов -на узлах кластера (см. рис. [[fig:subord-ppl]]). +на узлах кластера (см. рис.\nbsp{}[[fig:subord-ppl]]). По принципу работу механизм управляющих объектов и конвейеров напоминает механизм работы процедур и стеков вызовов, с тем лишь преимуществом, что методы объектов вызываются асинхронно и параллельно друг другу (насколько это позволяет -логика программы). Поля управляющего объекта --- это локальные переменные стека, -метод ~act~ --- это последовательность процессорных инструкций перед вложенным -вызовом процедуры, а метод ~react~ --- это последовательность инструкций после -вложенного вызова. Создание и отправка на конвейер подчиненного объекта --- это +логика программы). Поля управляющего объекта\nbsp{}--- это локальные переменные стека, +метод ~act~\nbsp{}--- это последовательность процессорных инструкций перед вложенным +вызовом процедуры, а метод ~react~\nbsp{}--- это последовательность инструкций после +вложенного вызова. Создание и отправка на конвейер подчиненного объекта\nbsp{}--- это вложенный вызов процедуры. Наличие двух методов обуславливается асинхронностью вложенных вызовов и помогает заменить активное ожидание завершения подчиненных объектов пассивным при помощи конвейеров. Конвейеры, в свою очередь, позволяют @@ -2167,13 +2151,13 @@ graph G { планировщик возвращает поток управления в родительский управляющий объект каждый раз когда какой-либо его дочерний объект завершает выполнение. Такое взаимодействие превращает сопрограмму в некоторого рода обработчик событий, в - котором событием является дочерний объект, а обработчиком --- родительский. + котором событием является дочерний объект, а обработчиком\nbsp{}--- родительский. - Сопрограмма может взаимодействовать с произвольным количеством управляющих объектов, адреса которых известны; взаимодействие с объектами, осуществляемое вразрез с иерархией сильно усложняет поток управления и стек вызовов сопрограмм теряет древовидную структуру. Только логика программы может гарантировать существование в памяти машины двух взаимодействующих объектов. - Один из способов обеспечения такой гарантии --- взаимодействие между + Один из способов обеспечения такой гарантии\nbsp{}--- взаимодействие между вложенными сопрограммами, вызванными из одной родительской сопрограммы. Поскольку такого рода взаимодействие можно осуществить в рамках иерархии через родительскую сопрограмму, его можно считать оптимизацией, позволяющей @@ -2234,8 +2218,7 @@ graph G { Для того чтобы учесть неоднородность частей, на которые разбиваются входные данные, и неоднородность выполняемых задач, необходимо предсказать время -выполнения каждой из задач. Соответствующее исследование сделано в -cite:degtyarev2016balance, поскольку реализация модели АРСС включает в себя, в +выполнения каждой из задач. Соответствующее исследование сделано в\nbsp{}cite:degtyarev2016balance, поскольку реализация модели АРСС включает в себя, в основном, однородные задачи. Таким образом, распределение нагрузки осуществляется в два этапа: на первом @@ -2253,10 +2236,10 @@ cite:degtyarev2016balance, поскольку реализация модели многопроцессорной машины достаточно для создания типовых реализаций морского волнения. Также использование видеокарт в качестве векторных ускорителей эффективно только в случае расчета давлений, в то время как генерация волновой -поверхности выполняется быстрее на скалярном процессоре cite:degtyarev2011effi. +поверхности выполняется быстрее на скалярном процессоре\nbsp{}cite:degtyarev2011effi. Создание программной реализации происходило в два этапа: на первом этапе был -создан и отлажен прототип в программной среде Mathematica cite:mathematica10, а +создан и отлажен прототип в программной среде Mathematica\nbsp{}cite:mathematica10, а на втором этапе логика программы была переписана на более низкоуровневом языке C++, и для получения эффективно работающего параллельного кода были проведены эксперименты с рядом библиотек. С помощью этих библиотек были реализованы @@ -2267,10 +2250,9 @@ C++, и для получения эффективно работающего п использование видеокарт неэффективно при генерации волновой поверхности (см. [[tab:autoreg-performance]]), что обусловлено сравнительно небольшим количеством арифметических операций по отношению к количеству операций с памятью устройства, -а также отсутствием трансцендентных функций в реализации алгоритма -cite:degtyarev2011effi. Во-вторых, для генерации одной реализации взволнованной +а также отсутствием трансцендентных функций в реализации алгоритма\nbsp{}cite:degtyarev2011effi. Во-вторых, для генерации одной реализации взволнованной морского поверхности одной многопроцессорной машины достаточно для эффективного -и быстрого решения задачи (см. [[fig:autoreg-performance]]). По результатам +и быстрого решения задачи (см. рис.\nbsp{}[[fig:autoreg-performance]]). По результатам тестирования стандарт OpenMP был выбран в качестве основного, как наиболее эффективный и наиболее подходящий для расчетов на многопроцессорной системе. @@ -2312,10 +2294,10 @@ cite:degtyarev2011effi. Во-вторых, для генерации одной была показана тестированием их разработчиками. В качестве библиотеки для матричных операций (расчета коэффициентов авторегрессионной модели) была выбрана GotoBLAS и основанная на ней LAPACK, для непрерывной аппроксимации поля волновых -чисел использовалась библиотека CGAL cite:fabri2009cgal и для статистической +чисел использовалась библиотека CGAL\nbsp{}cite:fabri2009cgal и для статистической проверки интегральных характеристик реализации взволнованной поверхности -использовалась библиотека GSL cite:gsl2008scientific. В случае GotoBLAS -эффективность библиотеки показана в работах cite:goto2008high,goto2008anatomy, +использовалась библиотека GSL\nbsp{}cite:gsl2008scientific. В случае GotoBLAS +эффективность библиотеки показана в работах\nbsp{}cite:goto2008high,goto2008anatomy, для других библиотек эффективность не является важной, и они были выбраны, исходя из удобства их использования. @@ -2324,12 +2306,12 @@ GotoBLAS и основанная на ней LAPACK, для непрерывно #+attr_latex: :booktabs t :align lp{0.6\linewidth} | Library | What it is used for | |--------------------------------------------------------+----------------------------------| -| DCMT cite:matsumoto1998dynamic | параллельный ГПСЧ | -| Blitz cite:veldhuizen1997will,veldhuizen2000techniques | многомерные массивы | -| GSL cite:gsl2008scientific | вычисление ФПР, ФР, БПФ | +| DCMT\nbsp{}cite:matsumoto1998dynamic | параллельный ГПСЧ | +| Blitz\nbsp{}cite:veldhuizen1997will,veldhuizen2000techniques | многомерные массивы | +| GSL\nbsp{}cite:gsl2008scientific | вычисление ФПР, ФР, БПФ | | | проверка стационарности процесса | -| LAPACK, GotoBLAS cite:goto2008high,goto2008anatomy | определение коэффициентов АР | -| GL, GLUT cite:kilgard1996opengl | трехмерная визуализация | +| LAPACK, GotoBLAS\nbsp{}cite:goto2008high,goto2008anatomy | определение коэффициентов АР | +| GL, GLUT\nbsp{}cite:kilgard1996opengl | трехмерная визуализация | **** Производительность алгоритма распределения нагрузки. Программная реализация генерации взволнованной поверхности сбалансирована с @@ -2370,9 +2352,9 @@ GotoBLAS и основанная на ней LAPACK, для непрерывно В эксперименте алгоритм распределения нагрузки показал большую эффективность по сравнению с реализацией без него. Чем больше размер генерируемой поверхности, -тем больше разрыв в производительности (рис. [[fig:factory-performance]]), что +тем больше разрыв в производительности (рис.\nbsp{}[[fig:factory-performance]]), что является следствием наложения вычислительной фазы и фазы вывода данных друг на -друга (рис. [[fig:factory-overlap]]). В реализации OpenMP фаза вывода данных +друга (рис.\nbsp{}[[fig:factory-overlap]]). В реализации OpenMP фаза вывода данных начинается только тогда, когда заканчивается вычислительная фаза, в то время как использование алгоритма распределения нагрузки приводит почти к одновременному завершению обеих фаз. Таким образом, /выполнение параллельных изнутри, @@ -2437,8 +2419,7 @@ arma.plot_factory_vs_openmp_overlap( требует наличия простаивающих резервных узлов на случай отказа главного узла. Алгоритмы выбора лидера (которые иногда называют алгоритмами /распределенного -консенсуса/) являются частными случаями волновых алгоритмов. В -cite:tel2000introduction Тель определяет их как алгоритмы, в которых событие +консенсуса/) являются частными случаями волновых алгоритмов. В\nbsp{}cite:tel2000introduction Тель определяет их как алгоритмы, в которых событие завершения программы предваряется хотя бы одним каким-либо другим событием, происходящем в /каждом/ параллельном процессе. Волновые алгоритмы не определены для анонимных сетей, т.е. они работают только с теми параллельными процессами, @@ -2470,7 +2451,7 @@ cite:tel2000introduction Тель определяет их как алгори IP-адресов подсети. Следующие ключевые особенности отличают наш подход от некоторых предложенных -ранее подходов cite:brunekreef1996design,aguilera2001stable,romano2014design. +ранее подходов\nbsp{}cite:brunekreef1996design,aguilera2001stable,romano2014design. - *Многоуровневая иерархия.* Количество руководящих узлов в сети зависит от значения ветвления. Если оно меньше количества IP-адресов в подсети, то в кластере будет несколько руководящих узлов. Если оно больше или равно @@ -2543,7 +2524,7 @@ IP-адреса: замена отображения IP-адресов на чт \forall f \colon \mathcal{N} \rightarrow \mathcal{R}^n \Rightarrow (f(n_1) < f(n_2) \Leftrightarrow \neg (f(n_1) \geq f(n_2))), \end{equation*} -где \(f\) --- отображение узла на его ранг, а \(<\) --- оператор, определяющий +где \(f\)\nbsp{}--- отображение узла на его ранг, а \(<\)\nbsp{}--- оператор, определяющий отношение строго порядка на множестве \(\mathcal{R}^n\). Функция \(f\) присваивает узлу порядковый номер, а оператор \(<\) делает этот номер уникальным. @@ -2570,7 +2551,7 @@ IP-адреса: замена отображения IP-адресов на чт l \geq 0, \quad o \geq 0 \end{equation*} -где \(n\) --- позиция IP-адреса узла в диапазоне IP-адресов подсети и \(p\) --- +где \(n\)\nbsp{}--- позиция IP-адреса узла в диапазоне IP-адресов подсети и \(p\)\nbsp{}--- значение ветвления (максимальное количество подчиненных, которых может иметь узел). Руководитель узла на уровне \(l\) с отступом \(o\) имеет уровень \(l-1\) и отступ \(\lfloor{o/p}\rfloor\). Расстояние между любыми двумя узлами в иерархии, @@ -2608,8 +2589,7 @@ IP-адреса: замена отображения IP-адресов на чт Платформа, на которой осуществлялось тестирование, представляла собой несколько многоядерных узлов, поверх которых с помощью пространств имен Linux развертывался виртуальный кластер из заданного количества узлов. Похожий подход -используется в -cite:lantz2010network,handigol2012reproducible,heller2013reproducible, где +используется в\nbsp{}cite:lantz2010network,handigol2012reproducible,heller2013reproducible, где авторы воспроизводят разнообразные практические эксперименты на виртуальных кластерах и сопоставляют результаты с физическими. Преимущество данного подхода заключается в возможности проведения экспериментов на больших виртуальных @@ -2627,8 +2607,9 @@ cite:lantz2010network,handigol2012reproducible,heller2013reproducible, где использованы дополнительные физические узлы, на каждом из которых запускалось по 100 виртуальных. Эксперимент показал, что обнаружение 100--400 узлами друг друга занимает в среднем 1,5 секунды, и это значение ненамного увеличивается с ростом -количества узлов (см. [[fig:bootstrap-local]]). Пример древовидной иерархии для 11 -узлов с ветвлением равным 2 представлен на [[fig:tree-hierarchy-11]]. +количества узлов (см.\nbsp{}рис.\nbsp{}[[fig:bootstrap-local]]). Пример древовидной +иерархии для 11 узлов с ветвлением равным 2 представлен на +рис.\nbsp{}[[fig:tree-hierarchy-11]]. #+name: fig:bootstrap-local #+caption: Зависимость времени объединения узлов в кластер от их количества. @@ -2688,7 +2669,7 @@ digraph { восстановления. Оптимизации работы контрольных точек восстановления посвящено большое -количество работ cite:egwutuoha2013survey, а альтернативным подходам +количество работ\nbsp{}cite:egwutuoha2013survey, а альтернативным подходам уделяется меньше внимания. Обычно высокопроизводительные приложения используют передачу сообщений для обмена данными между параллельными процессами и хранят свое текущее состояние в глобальной памяти, поэтому не существует способа @@ -2696,7 +2677,7 @@ digraph { памяти на диск. Обычно общее число процессов фиксировано и задается планировщиком, и в случае отказа перезапускаются сразу все процессы. Существуют некоторые обходные решения, которые позволяют перезапустить только часть -процессов cite:meyer2012radic, восстановив их на других узлах, однако это +процессов\nbsp{}cite:meyer2012radic, восстановив их на других узлах, однако это может привести к перегрузке, если на этих узлах уже запущены другие задачи. Теоретически, перезапуск процесса необязателен если задача может быть продолжена на меньшем количестве узлов, но библиотека передачи сообщений не @@ -2712,13 +2693,12 @@ digraph { случае нагрузка должна быть динамически перераспределена между оставшимися узлами. Несмотря на то что динамическое распределение нагрузки было реализовано в поверх библиотеки передачи сообщений в ряде -работ cite:bhandarkar2001adaptive,lusk2010more, оно никогда не применялось в -задаче обеспечения отказоустойчивости. В этом разделе исследуются методы -обеспечения отказоустойчивости при выходе из строя подчиненных и главных узлов -и показывается, как приемы объектно-ориентированного программирования могут -быть использованы для сохранения минимального состояния программы, необходимого -для ее перезапуска, в иерархии объектов, а не в глобальных и локальных -переменных. +работ\nbsp{}cite:bhandarkar2001adaptive,lusk2010more, оно никогда не применялось +в задаче обеспечения отказоустойчивости. В этом разделе исследуются методы +обеспечения отказоустойчивости при выходе из строя подчиненных и главных узлов и +показывается, как приемы объектно-ориентированного программирования могут быть +использованы для сохранения минимального состояния программы, необходимого для +ее перезапуска, в иерархии объектов, а не в глобальных и локальных переменных. **** Иерархия управляющих объектов. Для распределения нагрузки узлы кластера объединяются в древовидную иерархию @@ -2728,7 +2708,7 @@ digraph { делает систему симметричной и легкой в обслуживании: на каждом узле установлен один и тот же набор программного обеспечения, что позволяет заменить один узел другим при выходе из строя первого. Похожее архитектурное решение используется в -хранилищах типа ключ-значение cite:anderson2010couchdb,lakshman2010cassandra для +хранилищах типа ключ-значение\nbsp{}cite:anderson2010couchdb,lakshman2010cassandra для обеспечения отказоустойчивости, однако автору неизвестны планировщики задач, которые используют данный подход. @@ -2737,9 +2717,8 @@ digraph { а дополнительные узлы используются по необходимости. Такое решение позволяет использовать произвольное количество узлов для запуска задачи и динамически менять это количество во время ее выполнения. Похожее решение используется в -системах обработки больших объемов данных -cite:dean2008mapreduce,vavilapalli2013yarn --- пользователь, запускающий задачу -на кластере, не указывает количество узлов, фактические узлы --- это узлы, на +системах обработки больших объемов данных\nbsp{}cite:dean2008mapreduce,vavilapalli2013yarn\nbsp{}--- пользователь, запускающий задачу +на кластере, не указывает количество узлов, фактические узлы\nbsp{}--- это узлы, на которых расположены входные файлы. С математической точки зрения управляющий объект \(K\) может быть определен как @@ -2754,7 +2733,7 @@ cite:dean2008mapreduce,vavilapalli2013yarn --- пользователь, зап используется для остановки рекурсии, и передается в качестве аргумента главному управляющему объекту программы. Аргумент управляющего объекта интерпретируется следующим образом. -- Если объект является только что созданным объектом, то аргумент --- это его +- Если объект является только что созданным объектом, то аргумент\nbsp{}--- это его родительский объект. - В остальных случаях аргументом может являться любой объект (чаще всего дочерний по отношению к текущему). @@ -2769,7 +2748,7 @@ cite:dean2008mapreduce,vavilapalli2013yarn --- пользователь, зап пула объекты могут быть переданы на другие узлы кластера без явного указания в исходном коде программы. -Вычислительные объекты реализованы в виде замыканий (функторы в C++) --- +Вычислительные объекты реализованы в виде замыканий (функторы в C++)\nbsp{}--- объектов-функций, которые сохраняют в себе аргументы, ссылку на породивший их объект и данные из предметной области задачи. Данные обрабатываются либо при выполнении объекта, либо для параллельной обработки создаются подчиненные @@ -2778,9 +2757,8 @@ cite:dean2008mapreduce,vavilapalli2013yarn --- пользователь, зап **** Обработка выхода узлов из строя. Наиболее распространенная стратегия при выходе из строя подчиненного узла -является перезапуск выполнявшихся на нем объектов на рабочих узлах --- -стратегия, которой следует язык Erlang при перезапуске подчиненных процессов -cite:armstrong2003thesis. Для того что реализовать этот метод в рамках иерархии +является перезапуск выполнявшихся на нем объектов на рабочих узлах\nbsp{}--- +стратегия, которой следует язык Erlang при перезапуске подчиненных процессов\nbsp{}cite:armstrong2003thesis. Для того что реализовать этот метод в рамках иерархии управляющих объектов, узел-отправитель сохраняет каждый объект, передаваемый на другие узлы кластера, и в случае отказа произвольного количества узлов, на которые были переданы объекты, их копии перераспределяются между оставшимися @@ -2933,15 +2911,15 @@ TODO translate сбоев, но с количеством узлов минус один. Это отношение получается из того же самого эксперимента и представлено на рис.\nbsp{}[[fig:slowdown]]. Разница в производительности в случае выхода из строя руководящего и подчиненного узлов -находится в пределах 5%, а в случае выхода из строя резервного узла --- в +находится в пределах 5%, а в случае выхода из строя резервного узла\nbsp{}--- в пределах 50% для количества узлов меньше 6[fn::Измерение разницы для большего количества узлов не имеет смысла, поскольку программа завершается еще до -наступления сбоя.]. Увеличение времени выполнения на 50% --- это больше, чем -\(1/3\) времени работы программы, после которого происходит сбой, однако отказ -резервного узла требует некоторого времени, чтобы быть обнаруженным другими -узлами: сбой узла обнаруживается только тогда, когда подчиненный объект, имеющий -копию главного объекта, завершает свое выполнение и пытается вернуться на -исходный узел к родителю. Мгновенное обнаружение сбоя узла требует внезапной +наступления сбоя.]. Увеличение времени выполнения на 50%\nbsp{}--- это больше, +чем \(1/3\) времени работы программы, после которого происходит сбой, однако +отказ резервного узла требует некоторого времени, чтобы быть обнаруженным +другими узлами: сбой узла обнаруживается только тогда, когда подчиненный объект, +имеющий копию главного объекта, завершает свое выполнение и пытается вернуться +на исходный узел к родителю. Мгновенное обнаружение сбоя узла требует внезапной остановки выполнения объектов, что может быть неприменимо для программ со сложной логикой. @@ -3077,13 +3055,13 @@ TODO translate * Благодарности Графики в этой работе были подготовлены с помощью языка для статистических -вычислений R cite:rlang2016,Sarkar2008lattice и программного обеспечения -Graphviz cite:Gansner00anopen. Документ был подготовлен с использованием -Org-mode cite:Schulte2011org2,Schulte2011org1,Dominik2010org для GNU Emacs, -предоставляющего вычислительное окружение для воспроизводимых исследований. Это -означает, что все графики можно воспроизвести и соответствующие утверждения -проверить, скопировав репозиторий диссертации[fn:repo], установив Emacs и -экспортировав документ. +вычислений R\nbsp{}cite:rlang2016,Sarkar2008lattice и программного обеспечения +Graphviz\nbsp{}cite:Gansner00anopen. Документ был подготовлен с использованием +Org-mode\nbsp{}cite:Schulte2011org2,Schulte2011org1,Dominik2010org для GNU +Emacs, предоставляющего вычислительное окружение для воспроизводимых +исследований. Это означает, что все графики можно воспроизвести и +соответствующие утверждения проверить, скопировав репозиторий +диссертации[fn:repo], установив Emacs и экспортировав документ. [fn:repo] [[https://github.com/igankevich/arma-thesis]] diff --git a/phd-diss.org b/phd-diss.org @@ -260,11 +260,10 @@ Software programmes, which simulates vessel behaviour in sea waves, are widely used to model ship motion, estimate impact of external forces on floating platform or other marine object, and estimate capsize probability under given weather conditions; however, to model ocean waves most of the simulation codes -use linear wave theory -cite:shin2003nonlinear,van2007forensic,kat2001prediction,van2002development, in +use linear wave theory\nbsp{}cite:shin2003nonlinear,van2007forensic,kat2001prediction,van2002development, in the framework of which it is difficult to reproduce certain peculiarities of wind wave climate. Among them are transition between normal and storm weather, -and sea composed of multiple wave systems --- both wind waves and swell --- +and sea composed of multiple wave systems\nbsp{}--- both wind waves and swell\nbsp{}--- heading from multiple directions. Another shortcoming of linear wave theory is an assumption, that wave amplitude is small compared to wave length. This makes calculations imprecise when modelling ship motion in irregular waves, for which @@ -308,7 +307,7 @@ required mathematical apparatus. and can be generated as steep as it is possible with real ocean wave ACF. 3. Period of the realisation equals the period of PRNG, so generation time grows linearly with the realisation size. -4. White noise --- the only probabilistic term in ARMA process --- has +4. White noise\nbsp{}--- the only probabilistic term in ARMA process\nbsp{}--- has Gaussian distribution; so, convergence rate is not probabilistic. **** Goals and objectives. @@ -368,8 +367,7 @@ software. **** Methodology and research methods. Software implementation of ARMA model and pressure field formula was created -incrementally: a prototype written in high-level engineering language -cite:mathematica10,octave2015 was rewritten in lower level language (C++). +incrementally: a prototype written in high-level engineering language\nbsp{}cite:mathematica10,octave2015 was rewritten in lower level language (C++). Implementation of the same algorithm and formulae in languages of varying levels (which involves usage of different abstractions and language primitives) allows correcting errors, which would left unnoticed otherwise. Wavy surface, @@ -393,7 +391,7 @@ Mathematica language, where resulting formulae are verified by built-in graphical means. ARMA model and pressure field formula were incorporated into Large Amplitude -Motion Programme (LAMP) --- an ship motion simulation software programme --- +Motion Programme (LAMP)\nbsp{}--- an ship motion simulation software programme\nbsp{}--- where they were compared to previously used LH model. Preliminary numerical experiments showed higher computational efficiency of ARMA model. @@ -412,16 +410,16 @@ wave theory. **** Pressure field formula. The problem of finding pressure field under wavy sea surface represents inverse problem of hydrodynamics for incompressible inviscid fluid. System of equations -for it in general case is written as cite:kochin1966theoretical +for it in general case is written as\nbsp{}cite:kochin1966theoretical \begin{align} & \nabla^2\phi = 0,\nonumber\\ & \phi_t+\frac{1}{2} |\vec{\upsilon}|^2 + g\zeta=-\frac{p}{\rho}, & \text{на }z=\zeta(x,y,t),\label{eq:problem}\\ & D\zeta = \nabla \phi \cdot \vec{n}, & \text{на }z=\zeta(x,y,t),\nonumber \end{align} -where \(\phi\) --- velocity potential, \(\zeta\) --- elevation (\(z\) coordinate) -of wavy surface, \(p\) --- wave pressure, \(\rho\) --- fluid density, -\(\vec{\upsilon}=(\phi_x,\phi_y,\phi_z)\) --- velocity vector, \(g\) --- -acceleration of gravity, and \(D\) --- substantial (Lagrange) derivative. The +where \(\phi\)\nbsp{}--- velocity potential, \(\zeta\)\nbsp{}--- elevation (\(z\) coordinate) +of wavy surface, \(p\)\nbsp{}--- wave pressure, \(\rho\)\nbsp{}--- fluid density, +\(\vec{\upsilon}=(\phi_x,\phi_y,\phi_z)\)\nbsp{}--- velocity vector, \(g\)\nbsp{}--- +acceleration of gravity, and \(D\)\nbsp{}--- substantial (Lagrange) derivative. The first equation is called continuity (Laplace) equation, the second one is the conservation of momentum law (the so called dynamic boundary condition); the third one is kinematic boundary condition for free wavy surface, which states @@ -434,7 +432,7 @@ for \(\phi\). In this formulation dynamic boundary condition becomes explicit formula to determine pressure field using velocity potential derivatives obtained from the remaining equations. So, from mathematical point of view inverse problem of hydrodynamics reduces to Laplace equation with mixed boundary -condition --- Robin problem. +condition\nbsp{}--- Robin problem. * Related work ** Ocean wave models analysis @@ -446,9 +444,8 @@ generation, instead of wavy surface generation. *** Longuet---Higgins model The simplest model, formula of which is derived in the framework of linear wave -theory, is Longuet---Higgins (LH) model cite:longuet1957statistical. In-depth -comparative analysis of this model and ARMA model is done in -cite:degtyarev2011modelling,boukhanovsky1997thesis. +theory, is Longuet---Higgins (LH) model\nbsp{}cite:longuet1957statistical. In-depth +comparative analysis of this model and ARMA model is done in\nbsp{}cite:degtyarev2011modelling,boukhanovsky1997thesis. LH model represents ocean wavy surface as a superposition of sine waves with random amplitudes \(c_n\) and phases \(\epsilon_n\), continuously @@ -466,7 +463,7 @@ Gaussian process defined by \begin{equation*} 2E_\zeta(u,v)\, du\, dv = \sum\limits_n c_n^2, \end{equation*} -where \(E_\zeta(u,v)\) --- two-dimensional wave energy spectral density. +where \(E_\zeta(u,v)\)\nbsp{}--- two-dimensional wave energy spectral density. Coefficients \(c_n\) are derived from wave energy spectrum \(S(\omega)\) via \begin{equation*} c_n = \sqrt{ \textstyle\int\limits_{\omega_n}^{\omega_{n+1}} S(\omega) d\omega}. @@ -481,8 +478,8 @@ appear in practice. amplitudes and phases has normal distribution, no matter what spectrum is used as the model input. Using lower number of coefficients may solve the problem, but also make realisation period smaller. So, using LH model to - simulate waves with non-Gaussian distribution of elevation --- a distribution - which real ocean waves have cite:huang1980experimental,рожков1996теория --- + simulate waves with non-Gaussian distribution of elevation\nbsp{}--- a distribution + which real ocean waves have\nbsp{}cite:huang1980experimental,рожков1996теория\nbsp{}--- is impractical. 2. From computational point of view, the deficiency of the model is non-linear increase of wavy surface generation time with the increase of realisation @@ -496,21 +493,20 @@ appear in practice. low due to randomness of phases \(\epsilon_n\). - It is difficult to generalise LH model for non-Gaussian processes as it involves incorporating non-linear terms in ([[eq:longuet-higgins]]) for which - there is no known formula to determine coefficients - cite:рожков1990вероятностные. + there is no known formula to determine coefficients\nbsp{}cite:рожков1990вероятностные. To summarise, LH model is linear, computationally inefficient for long-time simulations, and difficult to use as a base for more advanced models. *** ARMA model -In cite:spanos1982arma ARMA model is used to generate time series spectrum of +In\nbsp{}cite:spanos1982arma ARMA model is used to generate time series spectrum of which is compatible with Pierson---Moskowitz (PM) approximation of ocean wave spectrum. The authors carry out experiments for one-dimensional AR, MA and ARMA models. They mention excellent agreement between target and initial spectra and higher performance of ARMA model compared to models based on summing large number of harmonic components with random phases. The also mention that in order to reach agreement between target and initial spectrum MA model require lesser -number of coefficients than AR model. In cite:spanos1996efficient the authors +number of coefficients than AR model. In\nbsp{}cite:spanos1996efficient the authors generalise ARMA model coefficients determination formulae for multi-variate (vector) case. @@ -533,7 +529,7 @@ opportunity to visualise output of the programme that allowed to ensure that generated surface is compatible with real ocean surface, and is not abstract multi-dimensional stochastic process that is real only statistically. -In cite:fusco2010short AR model is used to predict swell waves to control +In\nbsp{}cite:fusco2010short AR model is used to predict swell waves to control wave-energy converters (WEC) in real-time. In order to make WEC more efficient its internal oscillator frequency should match the one of ocean waves. The authors treat wave elevation as time series and compare performance of AR model, @@ -544,7 +540,7 @@ process to ocean wave modelling. ** Pressure field determination formulae *** Small amplitude waves theory -In cite:stab2012,детярев1998моделирование,degtyarev1997analysis the authors +In\nbsp{}cite:stab2012,детярев1998моделирование,degtyarev1997analysis the authors propose a solution for inverse problem of hydrodynamics of potential flow in the framework of small-amplitude wave theory (under assumption that wave length is much larger than height: \(\lambda \gg h\)). In that case inverse problem is @@ -625,9 +621,8 @@ distributed random impulses. The governing equation for 3-D ARMA process is , \label{eq:arma-process} \end{equation} -where \(\zeta\) --- wave elevation, \(\Phi\) --- AR process coefficients, \(\Theta\) ---- MA process coefficients, \(\epsilon\) --- white noise with Gaussian -distribution, \(\vec N\) --- AR process order, \(\vec M\) --- MA process order, and +where \(\zeta\)\nbsp{}--- wave elevation, \(\Phi\)\nbsp{}--- AR process coefficients, \(\Theta\)\nbsp{}--- MA process coefficients, \(\epsilon\)\nbsp{}--- white noise with Gaussian +distribution, \(\vec N\)\nbsp{}--- AR process order, \(\vec M\)\nbsp{}--- MA process order, and \(\Phi_{\vec{0}}\equiv{0}\), \(\Theta_{\vec{0}}\equiv{0}\). Here arrows denote multi-component indices with a component for each dimension. In general, any scalar quantity can be a component (temperature, salinity, concentration of some @@ -667,7 +662,7 @@ Generic form of YW equations is 0, \quad \text{if } \vec{k}\neq0, \end{cases} \end{equation} -where \(\gamma\) --- ACF of process \(\zeta\), \(\Var{\epsilon}\) --- white noise +where \(\gamma\)\nbsp{}--- ACF of process \(\zeta\), \(\Var{\epsilon}\)\nbsp{}--- white noise variance. Matrix form of three-dimensional YW equations, which is used in the present work, is \begin{equation*} @@ -778,7 +773,7 @@ Here coefficients \(\Theta\) are calculated from back to front: from \Theta_{\vec j}^2 }. \end{equation*} -Authors of cite:box1976time suggest using Newton---Raphson method to solve this +Authors of\nbsp{}cite:box1976time suggest using Newton---Raphson method to solve this equation with higher precision, however, this method does not work in three dimensions. Using slower method does not have dramatic effect on the overall programme performance, because the number of coefficients is small and most of @@ -793,7 +788,7 @@ Generally speaking, ARMA process is obtained by plugging MA generated wavy surface as random impulse to AR process, however, in order to get the process with desired ACF one should re-compute AR coefficients before plugging. There are several approaches to "mix" AR and MA processes. -- The approach proposed in cite:box1976time which involves dividing ACF into MA +- The approach proposed in\nbsp{}cite:box1976time which involves dividing ACF into MA and AR part along each dimension is not applicable here, because in three dimensions such division is not possible: there always be parts of the ACF that are not taken into account by AR and MA process. @@ -817,13 +812,13 @@ use AR process for standing waves and MA process for progressive waves. The other problem is inability to automatically determine optimal number of coefficients for three-dimensional AR and MA processes. For one-dimensional -processes this can be achieved via iterative methods cite:box1976time, but they +processes this can be achieved via iterative methods\nbsp{}cite:box1976time, but they diverge in three-dimensional case. The final problem, which is discussed in [[#sec:how-to-mix-ARMA]], is inability to "mix" AR and MA process in three dimensions. -In practice some statements made for AR and MA processes in cite:box1976time +In practice some statements made for AR and MA processes in\nbsp{}cite:box1976time should be flipped for three-dimensional case. For example, the authors say that ACF of MA process cuts at \(q\) and ACF of AR process decays to nought infinitely, but in practice making ACF of 3-dimensional MA process not decay results in it @@ -844,18 +839,17 @@ future research. ** Modelling non-linearity of ocean waves ARMA model allows modelling asymmetry of wave elevation distribution, i.e. generate ocean waves, distribution of z-coordinate of which has non-nought -kurtosis and asymmetry. Such distribution is inherent to real ocean waves -cite:longuet1963nonlinear. +kurtosis and asymmetry. Such distribution is inherent to real ocean waves\nbsp{}cite:longuet1963nonlinear. Wave asymmetry is modelled by non-linear inertia-less transform (NIT) of stochastic process, however, transforming resulting wavy surface means transforming initial ACF. In order to alleviate this, ACF must be preliminary -transformed as shown in cite:boukhanovsky1997thesis. +transformed as shown in\nbsp{}cite:boukhanovsky1997thesis. **** Wavy surface transformation. Explicit formula \(z=f(y)\) that transforms wavy surface to desired one-dimensional distribution \(F(z)\) is the solution of non-linear transcendental -equation \(F(z)=\Phi(y)\), where \(\Phi(y)\) --- one-dimensional Gaussian +equation \(F(z)=\Phi(y)\), where \(\Phi(y)\)\nbsp{}--- one-dimensional Gaussian distribution. Since distribution of wave elevation is often given by some approximation based on field data, this equation is solved numerically with respect to \(z_k\) in each grid point \(y_k|_{k=0}^N\) of generated wavy surface. In @@ -886,7 +880,7 @@ where \int\limits_{0}^\infty f(y) H_m(y) \exp\left[ -\frac{y^2}{2} \right], \end{equation*} -\(H_m\) --- Hermite polynomial, and \(f(y)\) --- solution to equation +\(H_m\)\nbsp{}--- Hermite polynomial, and \(f(y)\)\nbsp{}--- solution to equation eqref:eq:distribution-transformation. Plugging polynomial approximation \(f(y)\approx\sum\limits_{i}d_{i}y^i\) and analytic formulae for Hermite polynomial yields @@ -908,7 +902,7 @@ fields become equal with desired accuracy \(\epsilon\): \frac{C_k^2}{k!} \right| \leq \epsilon. \end{equation*} -In cite:boukhanovsky1997thesis the author suggests using polynomial +In\nbsp{}cite:boukhanovsky1997thesis the author suggests using polynomial approximation \(f(y)\) also for wavy surface transformation, however, in practice ocean surface realisation often contains points, where z-coordinate is beyond the limits of the approximation, which makes solution wrong. In these points it @@ -941,7 +935,7 @@ comparable to that of FFT. For example, stationary elliptic PDE transforms to implicit numerical scheme which is solved by iterative method on each step of which a tridiagonal of five-diagonal system of algebraic equations is solved by Thomas algorithm. Asymptotic complexity of this approach is -\(\mathcal{O}({n}{m})\), where \(n\) --- number of wavy surface grid points, \(m\) --- +\(\mathcal{O}({n}{m})\), where \(n\)\nbsp{}--- number of wavy surface grid points, \(m\)\nbsp{}--- number of iterations. Despite their wide spread, iterative algorithms are inefficient on parallel computer architectures; in particular, their mapping to co-processors may involve copying data in and out of the co-processor in each @@ -998,7 +992,7 @@ transforms, we rewrite eqref:eq:guessed-sol-2d as a convolution: \ast \InverseFourierY{E(u)}{x}, \end{equation*} -where \(\Fun{z}\) --- a function, form of which is defined in section +where \(\Fun{z}\)\nbsp{}--- a function, form of which is defined in section [[#sec:compute-delta]] and which satisfies equation \(\FourierY{\Fun{z}}{u}=e^{2\pi{u}{z}}\). Plugging formula \(\phi\) into the boundary condition yields @@ -1052,7 +1046,7 @@ following section. **** Formula for finite depth fluid. On the sea bottom vertical fluid velocity component equals nought: \(\phi_z=0\) on -\(z=-h\), where \(h\) --- water depth. In this case equation \(v=-{i}{u}\), which came +\(z=-h\), where \(h\)\nbsp{}--- water depth. In this case equation \(v=-{i}{u}\), which came from Laplace equation, can not be neglected, hence the solution is sought in the following form: \begin{equation} @@ -1099,14 +1093,14 @@ previous section transformations yields final formula for \(\phi(x,z)\): } \label{eq:solution-2d-full} \end{equation} -where \(\FunSecond{z}\) --- a function, form of which is defined in section +where \(\FunSecond{z}\)\nbsp{}--- a function, form of which is defined in section [[#sec:compute-delta]] and which satisfies equation \(\FourierY{\FunSecond{z}}{u}=\Sinh{2\pi{u}{z}}\). **** Reducing to the formulae from linear wave theory. Check the validity of derived formulae by substituting \(\zeta(x,t)\) with known analytic formula for plain waves. Symbolic computation of Fourier transforms in -this section were performed in Mathematica cite:mathematica10. In the framework +this section were performed in Mathematica\nbsp{}cite:mathematica10. In the framework of linear wave theory assume that waves have small amplitude compared to their lengths, which allows us simplifying initial system of equations eqref:eq:problem-2d to @@ -1251,13 +1245,13 @@ For three-dimensional wave profile (2D in space and 1D in time) analytic formula is a polynomial of high order and is best obtained via symbolic computation programme. Then for practical usage it can be approximated by superposition of exponentially decaying cosines (which is how ACF of a stationary ARMA process -looks like cite:box1976time). +looks like\nbsp{}cite:box1976time). **** Empirical method of finding the ACF. However, for three-dimensional case there exists simpler empirical method which does not require sophisticated software to determine shape of the ACF. It is known that ACF represented by exponentially decaying cosines satisfies first -order Stokes' equations for gravity waves cite:boccotti1983wind. So, if the +order Stokes' equations for gravity waves\nbsp{}cite:boccotti1983wind. So, if the shape of the wave profile is the only concern in the simulation, then one can simply multiply it by a decaying exponent to get appropriate ACF. This ACF does not reflect other wave profile parameters, such as wave height and period, but @@ -1353,9 +1347,9 @@ function (PDF) of the surface elevation. This distribution is given by either polynomial approximation of /in situ/ data or analytic formula. **** Gram---Charlier series expansion. -In cite:huang1980experimental the authors experimentally show, that PDF of sea +In\nbsp{}cite:huang1980experimental the authors experimentally show, that PDF of sea surface elevation is distinguished from normal distribution by non-nought -kurtosis and skewness. In cite:рожков1996теория the authors show, that this type +kurtosis and skewness. In\nbsp{}cite:рожков1996теория the authors show, that this type of PDF expands in Gram---Charlier series: \begin{align} \label{eq:skew-normal-1} @@ -1379,12 +1373,12 @@ of PDF expands in Gram---Charlier series: +1 \right], \end{align} -where \(\phi(z)=\frac{1}{2}\mathrm{erf}(z/\sqrt{2})\), \(\gamma_1\) --- skewness, -\(\gamma_2\) --- kurtosis, \(f\) --- PDF, \(F\) --- cumulative distribution function -(CDF). According to cite:рожков1990вероятностные for ocean waves skewness is +where \(\phi(z)=\frac{1}{2}\mathrm{erf}(z/\sqrt{2})\), \(\gamma_1\)\nbsp{}--- skewness, +\(\gamma_2\)\nbsp{}--- kurtosis, \(f\)\nbsp{}--- PDF, \(F\)\nbsp{}--- cumulative distribution function +(CDF). According to\nbsp{}cite:рожков1990вероятностные for ocean waves skewness is selected from interval \(0.1\leq\gamma_1\leq{0.52}]\) and kurtosis from interval \(0.1\leq\gamma_2\leq{0.7}\). Family of probability density functions for -different parameters is shown in [[fig:skew-normal-1]]. +different parameters is shown in fig.\nbsp{}[[fig:skew-normal-1]]. #+name: fig:skew-normal-1 #+begin_src R :results output graphics :exports results :file build/skew-normal-1.pdf @@ -1426,12 +1420,12 @@ elevation by skew-normal distribution: f(z; \alpha) & = \frac{e^{-\frac{z^2}{2}}}{\sqrt{2 \pi }} \mathrm{erfc}\left[-\frac{\alpha z}{\sqrt{2}}\right], \end{align} -where \(T\) --- Owen \(T\)-function cite:owen1956tables. Using this formula it is -impossible to specify skewness and kurtosis separately --- both values are +where \(T\)\nbsp{}--- Owen \(T\)-function\nbsp{}cite:owen1956tables. Using this formula it is +impossible to specify skewness and kurtosis separately\nbsp{}--- both values are adjusted via \(\alpha\) parameter. The only advantage of the formula is its relative computational simplicity: this function is available in some programmes and mathematical libraries. Its graph for different values of \(\alpha\) is shown -in [[fig:skew-normal-2]]. +in fig.\nbsp{}[[fig:skew-normal-2]]. #+name: fig:skew-normal-2 #+begin_src R :results output graphics :exports results :file build/skew-normal-2.pdf @@ -1483,14 +1477,13 @@ transform ACF; relative error without interpolation is \(10^{-5}\). *** White noise generation In order to eliminate periodicity from generated wavy surface, it is imperative to use PRNG with sufficiently large period to generate white noise. Parallel -Mersenne Twister cite:matsumoto1998mersenne with a period of \(2^{19937}-1\) is +Mersenne Twister\nbsp{}cite:matsumoto1998mersenne with a period of \(2^{19937}-1\) is used as a generator in this work. It allows producing aperiodic ocean wavy surface realisations in any practical usage scenarios. There is no guarantee that multiple Mersenne Twisters executed in parallel threads with distinct initial states produce uncorrelated pseudo-random number -sequences, however, algorithm of dynamic creation of Mersenne Twisters -cite:matsumoto1998dynamic may be used to provide such guarantee. The essence of +sequences, however, algorithm of dynamic creation of Mersenne Twisters\nbsp{}cite:matsumoto1998dynamic may be used to provide such guarantee. The essence of the algorithm is to find matrices of initial generator states, that give maximally uncorrelated pseudo-random number sequences when Mersenne Twisters are executed in parallel with these initial states. Since finding such initial @@ -1501,7 +1494,7 @@ saved to a file, which is then read before starting white noise generation. *** Wavy surface generation In ARMA model value of wavy surface elevation at a particular point depends on previous in space and time points, as a result the so called /ramp-up interval/ -(see fig. [[fig:ramp-up-interval]]), in which realisation does not correspond to +(see fig.\nbsp{}[[fig:ramp-up-interval]]), in which realisation does not correspond to specified ACF, forms in the beginning of the realisation. There are several solutions to this problem which depend on the simulation context. @@ -1518,8 +1511,7 @@ LH model and generate the rest of the realisation with ARMA model. Algorithm of wavy surface generation is data-parallel: realisation is divided into equal parts each of which is generated independently, however, in the beginning of each realisation there is ramp-up interval. To eliminate it -/overlap-add/ method -cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms (a popular +/overlap-add/ method\nbsp{}cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms (a popular method in signal processing) is used. The essence of the method is to add another interval, size of which is equal to the ramp-up interval size, to the end of each part. Then wavy surface is generated in each point of each part @@ -1573,7 +1565,7 @@ which terms with \(\zeta\) are omitted. :CUSTOM_ID: sec:verification :END: -In cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis AR +In\nbsp{}cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis AR model the following items are verified experimentally: - probability distributions of different wave characteristics (wave heights, lengths, crests, periods, slopes, three-dimensionality), @@ -1583,7 +1575,7 @@ In this work both AR and MA model are verified by comparing probability distributions of different wave characteristics. *** Verification of wavy surface integral characteristics -In cite:рожков1990вероятностные the authors show that several ocean wave +In\nbsp{}cite:рожков1990вероятностные the authors show that several ocean wave characteristics (listed in table [[tab:weibull-shape]]) have Weibull distribution, and wavy surface elevation has Gaussian distribution. In order to verify that distributions corresponding to generated realisation are correct, @@ -1680,7 +1672,7 @@ of [[#sec:pressure-2d]], only finite depth formulae are compared. The experiment shows that velocity potential fields produced by formula eqref:eq:solution-2d-full for finite depth fluid and formula eqref:eq:solution-2d-linear from linear wave theory are qualitatively different -(fig. [[fig:potential-field-nonlinear]]). First, velocity potential contours have +(fig.\nbsp{}[[fig:potential-field-nonlinear]]). First, velocity potential contours have sinusoidal shape, which is different from oval shape described by linear wave theory. Second, velocity potential decays more rapidly than in linear wave theory as getting closer to the bottom, and the region where the majority of @@ -1706,8 +1698,7 @@ Integration in formula eqref:eq:solution-2d-full is done over wave numbers range extracted from the generated wavy surface. For small-amplitude waves both formulae show comparable results (the difference in the velocity is attributed to stochastic nature of AR model), whereas for large-amplitude waves stable -velocity field is produced only by formula eqref:eq:solution-2d-full (fig. -[[fig:velocity-field-2d]]). So, generic formula eqref:eq:solution-2d-full gives +velocity field is produced only by formula eqref:eq:solution-2d-full (fig.\nbsp{}[[fig:velocity-field-2d]]). So, generic formula eqref:eq:solution-2d-full gives satisfactory results without restriction on wave amplitudes. #+name: fig:velocity-field-2d @@ -1735,7 +1726,7 @@ Software implementation of ARMA model works as a computational pipeline, in which each joint applies some function to the output coming from the pipe of the previous joint. Joints are distributed across computer cluster nodes to enable function parallelism, and then data flowing through the joints is distributed -across processor cores to enable data parallelism. Figure [[fig:pipeline]] shows a +across processor cores to enable data parallelism. Figure\nbsp{}[[fig:pipeline]] shows a diagram of data processing pipeline in which rectangles with rounded corners denote joints, regular rectangles denote arrays of problem domain objects flowing from one joint to another, and arrows show flow direction. Some joints @@ -1856,8 +1847,7 @@ digraph { [[file:build/pipeline.pdf]] Object pipeline may be seen as an improvement of BSP (Bulk Synchronous Parallel) -model cite:valiant1990bridging, which is used in graph processing -cite:malewicz2010pregel,seo2010hama. Pipeline eliminates global synchronisation +model\nbsp{}cite:valiant1990bridging, which is used in graph processing\nbsp{}cite:malewicz2010pregel,seo2010hama. Pipeline eliminates global synchronisation (where it is possible) after each sequential computation step by doing data transfer between joints in parallel to computations, whereas in BSP model global synchronisation occurs after each step. @@ -1907,8 +1897,7 @@ When one has such dependency, it is trivial to determine which task should be responsible for re-executing a failed task on one of the survived nodes. To re-execute the task on the top of the hierarchy, a backup task is created and executed on a different node. There exists a number of systems that are capable -of executing directed acyclic graphs of tasks in parallel -cite:acun2014charmpp,islam2012oozie, but graphs are not suitable to infer +of executing directed acyclic graphs of tasks in parallel\nbsp{}cite:acun2014charmpp,islam2012oozie, but graphs are not suitable to infer principal-subordinate relationship between tasks, because a node in the graph may have multiple parent nodes. @@ -1919,7 +1908,7 @@ it transparently to a programmer. The implementation is divided into two layers: the lower layer consists of routines and classes for single node applications (with no network interactions), and the upper layer for applications that run on an arbitrary number of nodes. There are two kinds of tightly coupled entities in -the model --- /control flow objects/ (or /kernels/) and /pipelines/ --- which +the model\nbsp{}--- /control flow objects/ (or /kernels/) and /pipelines/\nbsp{}--- which are used together to compose a programme. Kernels implement control flow logic in theirs ~act~ and ~react~ methods and @@ -1943,7 +1932,7 @@ thread pool that processes kernels in accordance with rules outlined in the previous paragraph. A separate pipeline is used for each device: There are pipelines for parallel processing, schedule-based processing (periodic and delayed tasks), and a proxy pipeline for processing of kernels on other cluster -nodes (see figure [[fig:subord-ppl]]). +nodes (see fig.\nbsp{}[[fig:subord-ppl]]). In principle, kernels and pipelines machinery reflect the one of procedures and call stacks, with the advantage that kernel methods are called asynchronously @@ -2143,8 +2132,7 @@ work in parallel to each other, lowering overall system resources downtime compared to using all devices from a single thread. In order to take into account non-homogeneous input data parts or tasks, one may -predict execution time of each task. Relevant study is done in -cite:degtyarev2016balance since ARMA model implementation includes mostly +predict execution time of each task. Relevant study is done in\nbsp{}cite:degtyarev2016balance since ARMA model implementation includes mostly homogeneous tasks. So, load balancing is done in two stages: in the first stage the task wrapped in @@ -2181,12 +2169,12 @@ efficient one. #+attr_latex: :booktabs t :align lp{0.6\linewidth} | Library | What it is used for | |--------------------------------------------------------+---------------------------------| -| DCMT cite:matsumoto1998dynamic | parallel PRNG | -| Blitz cite:veldhuizen1997will,veldhuizen2000techniques | multidimensional arrays | -| GSL cite:gsl2008scientific | PDF, CDF, FFT computation | +| DCMT\nbsp{}cite:matsumoto1998dynamic | parallel PRNG | +| Blitz\nbsp{}cite:veldhuizen1997will,veldhuizen2000techniques | multidimensional arrays | +| GSL\nbsp{}cite:gsl2008scientific | PDF, CDF, FFT computation | | | checking process stationarity | -| LAPACK, GotoBLAS cite:goto2008high,goto2008anatomy | finding AR coefficients | -| GL, GLUT cite:kilgard1996opengl | three-dimensional visualisation | +| LAPACK, GotoBLAS\nbsp{}cite:goto2008high,goto2008anatomy | finding AR coefficients | +| GL, GLUT\nbsp{}cite:kilgard1996opengl | three-dimensional visualisation | **** Performance of load balancing algorithm. Software implementation of wavy surface generation is balanced in terms of the @@ -2225,9 +2213,8 @@ and CPU thread pool size was equal the number of physical processor cores. In the experiment load balancing algorithm showed higher performance than implementation without it. The more the size of the generated surface is the -more the gap in performance is (fig. [[fig:factory-performance]]) which is a result -of overlap of computation phase and data output phase (fig. -[[fig:factory-overlap]]). In OpenMP implementation data output phase begins only +more the gap in performance is (fig.\nbsp{}[[fig:factory-performance]]) which is a result +of overlap of computation phase and data output phase (fig.\nbsp{}[[fig:factory-overlap]]). In OpenMP implementation data output phase begins only when computation is over, whereas load balancing algorithm makes both phases end almost simultaneously. So, /pipelined execution of internally parallel sequential phases is more efficient than their sequential execution/, and this @@ -2288,8 +2275,7 @@ algorithm, this approach becomes more and more popular as it does not require spare reserved nodes to recover from principal node failure. Leader election algorithms (which sometimes referred to as /distributed -consensus/ algorithms are special cases of wave algorithms. In -cite:tel2000introduction Tel defines them as algorithms in which termination +consensus/ algorithms are special cases of wave algorithms. In\nbsp{}cite:tel2000introduction Tel defines them as algorithms in which termination event is preceded by at least one event occurring in /each/ parallel process. Wave algorithms are not defined for anonymous networks, that is they apply only to processes that can uniquely define themselves. However, the number of @@ -2317,8 +2303,7 @@ from one node to another, which makes the cluster unmanageable. The simplest such function is the position of an IP address in network IP address range. The following key features distinguish this approach with respect to some -existing proposals -cite:brunekreef1996design,aguilera2001stable,romano2014design. +existing proposals\nbsp{}cite:brunekreef1996design,aguilera2001stable,romano2014design. - *Multi-level hierarchy.* The number of principal nodes in a network depends on the fan-out value. If it is lesser than the number of IP-addresses in the network, then there are multiple principle nodes in the cluster. If it is @@ -2445,8 +2430,7 @@ principal only, and inefficient scan of all network by each node does not occur. **** Evaluation results. Test platform consisted of several multi-core nodes, on top of which virtual clusters with varying number of nodes were deployed using Linux network -namespaces. Similar approach is used in -cite:lantz2010network,handigol2012reproducible,heller2013reproducible where the +namespaces. Similar approach is used in\nbsp{}cite:lantz2010network,handigol2012reproducible,heller2013reproducible where the authors reproduce various real-world experiments using virtual clusters and compare results to physical ones. The advantage of it is that the tests can be performed on a large virtual cluster using relatively small number of physical @@ -2462,9 +2446,9 @@ running more than 100 virtual nodes on one physical node simultaneously warp the results, thus additional physical nodes, each of which run 100 virtual nodes, were used for the experiment. The experiment showed that discovery of 100--400 nodes each other takes 1.5 seconds on average, and the value increases only -slightly with increase in the number of nodes (see fig. [[fig:bootstrap-local]]). An -example of tree hierarchy for 11 nodes with fan-out 2 is shown in fig. -[[fig:tree-hierarchy-11]]. +slightly with increase in the number of nodes (see +fig.\nbsp{}[[fig:bootstrap-local]]). An example of tree hierarchy for 11 nodes with +fan-out 2 is shown in fig.\nbsp{}[[fig:tree-hierarchy-11]]. #+name: fig:bootstrap-local #+caption: Time to discover all nodes of the cluster in depending on number of nodes. @@ -2540,16 +2524,15 @@ as a cluster operating system overlay allowing to write distributed applications. **** Related work. -Dynamic role assignment is an emerging trend in design of distributed systems -cite:ostrovsky2015couchbase,divya2013elasticsearch,boyer2012glusterfs,anderson2010couchdb,lakshman2010cassandra, +Dynamic role assignment is an emerging trend in design of distributed systems\nbsp{}cite:ostrovsky2015couchbase,divya2013elasticsearch,boyer2012glusterfs,anderson2010couchdb,lakshman2010cassandra, however, it is still not used in big data and HPC job schedulers. For example, -in popular YARN job scheduler cite:vavilapalli2013yarn, which is used by Hadoop +in popular YARN job scheduler\nbsp{}cite:vavilapalli2013yarn, which is used by Hadoop and Spark big data analysis frameworks, principal and subordinate roles are static. Failure of a subordinate node is tolerated by restarting a part of a job on a healthy node, and failure of a principal node is tolerated by setting up -standby reserved server cite:murthy2011architecture. Both principal servers are +standby reserved server\nbsp{}cite:murthy2011architecture. Both principal servers are coordinated by Zookeeper service which itself uses dynamic role assignment to -ensure its fault-tolerance cite:okorafor2012zookeeper. So, the whole setup is +ensure its fault-tolerance\nbsp{}cite:okorafor2012zookeeper. So, the whole setup is complicated due to Hadoop scheduler lacking dynamic roles: if dynamic roles were available, Zookeeper would be redundant in this setup. Moreover, this setup does not guarantee continuous operation of principal node because standby server @@ -2567,8 +2550,7 @@ it is far from ideal solution where roles are completely decoupled from physical servers. Finally, the simplest principal node high availability is implemented in Virtual -Router Redundancy Protocol (VRRP) -cite:knight1998rfc2338,hinden2004virtual,nadas2010rfc5798. Although VRRP +Router Redundancy Protocol (VRRP)\nbsp{}cite:knight1998rfc2338,hinden2004virtual,nadas2010rfc5798. Although VRRP protocol does provide principal and backup node roles, which are dynamically assigned to available routers, this protocol works on top of the IPv4 and IPv6 protocols and is designed to be used by routers and reverse proxy servers. Such @@ -2580,8 +2562,7 @@ daemon\nbsp{}cite:cassen2002keepalived. In contrast to web servers and HPC and big data job schedulers, some distributed key-value stores and parallel file systems have symmetric architecture, where principal and subordinate roles are assigned dynamically, so that any node can act as a -principal when the current principal node fails -cite:ostrovsky2015couchbase,divya2013elasticsearch,boyer2012glusterfs,anderson2010couchdb,lakshman2010cassandra. +principal when the current principal node fails\nbsp{}cite:ostrovsky2015couchbase,divya2013elasticsearch,boyer2012glusterfs,anderson2010couchdb,lakshman2010cassandra. This design decision simplifies management and interaction with a distributed system. From system administrator point of view it is much simpler to install the same software stack on each node than to manually configure principal and subordinate @@ -2610,18 +2591,17 @@ neighbours: when one runs the kernel on the subordinate node, the principal node also receive some of its subordinate kernels. This makes the system symmetrical and easy to maintain: each node have the same set of software that allows replacing one node with another in case of failure of the former. Similar -architectural solution used in key-value stores -cite:anderson2010couchdb,lakshman2010cassandra to provide fault tolerance, but +architectural solution used in key-value stores\nbsp{}cite:anderson2010couchdb,lakshman2010cassandra to provide fault tolerance, but author does not know any task schedulers that use this approach. Unlike ~main~ function in programmes based on message passing library, the first (the main) kernel is initially run only on one node, and remote nodes are used on-demand. This design choice allows having arbitrary number of nodes throughout execution of a programme, and use more nodes for highly parallel parts of the -code. Similar choice is made in the design of big data frameworks -cite:dean2008mapreduce,vavilapalli2013yarn --- a user submitting a job does not -specify the number of hosts to run its job on, and actual hosts are the hosts -where input files are located. +code. Similar choice is made in the design of big data +frameworks\nbsp{}cite:dean2008mapreduce,vavilapalli2013yarn\nbsp{}--- a user +submitting a job does not specify the number of hosts to run its job on, and +actual hosts are the hosts where input files are located. From mathematical point of view kernel \(K\) can be described as a vector-valued functional which recursively maps a kernel to \(n\)-component vector of kernels: @@ -2646,7 +2626,7 @@ unspecified order, several concurrent threads retrieve kernels from the pool and may send the remaining kernels to neighbouring cluster nodes if the pool overflows. -Kernels are implemented as closures (functors in C++) --- function objects +Kernels are implemented as closures (functors in C++)\nbsp{}--- function objects containing all their arguments, a reference to parent kernel and application domain data. The data is either processed upon kernel call, or subordinate kernels are created to process it in parallel. When the processing is complete a @@ -2655,8 +2635,8 @@ collect the resulting data from it. **** Handling nodes failures. Basic strategy to overcome a failure of a subordinate node is to restart -corresponding kernels on a healthy node --- a strategy employed by Erlang -language to restart failed subordinate processes cite:armstrong2003thesis. In +corresponding kernels on a healthy node\nbsp{}--- a strategy employed by Erlang +language to restart failed subordinate processes\nbsp{}cite:armstrong2003thesis. In order to implement this method in the framework of kernel hierarchy, sender node saves every kernel that is sent to remote cluster nodes, and in an event of a failure of any number of nodes, where kernels were sent, their copies are @@ -2923,12 +2903,12 @@ failures\nbsp{}cite:fekete1993impossibility. * Conclusion * Acknowledgements The graphs in this work were prepared using R language for statistical -computing cite:rlang2016,Sarkar2008lattice and Graphviz software -cite:Gansner00anopen. The manuscript was prepared using Org-mode -cite:Schulte2011org2,Schulte2011org1,Dominik2010org for GNU Emacs which provides -computing environment for reproducible research. This means that all graphs can -be reproduced and corresponding statements verified by cloning thesis -repository[fn:repo], installing Emacs and exporting the document. +computing\nbsp{}cite:rlang2016,Sarkar2008lattice and Graphviz +software\nbsp{}cite:Gansner00anopen. The manuscript was prepared using +Org-mode\nbsp{}cite:Schulte2011org2,Schulte2011org1,Dominik2010org for GNU Emacs +which provides computing environment for reproducible research. This means that +all graphs can be reproduced and corresponding statements verified by cloning +thesis repository[fn:repo], installing Emacs and exporting the document. [fn:repo] [[https://github.com/igankevich/arma-thesis]]