Об автоматическом дифференцировании (АД) на Хабре уже писалось здесь и здесь. В данной статье предлагается реализация АД для Delphi (протестировано в Embarcadero XE2, XE6) вместе с удобными классами методов Ньютона для решения нелинейных уравнений f(x) = 0 и систем F(X) = 0. Любые ссылки на готовые аналогичные библиотеки приветствуются, сам же я подобного не нашел, не считая отличного решателя СЛАУ с разряженной матрицей (см. под катом).
Давайте в самом начале отметим для себя, что выбор Delphi обусловлен legacy-кодом, тем не менее на C++ задачу можно решать следующим образом. Во-первых, для методов Ньютоновского (базовый метод Ньютона, метод Бройдена) типа имеются Numerical Recipes in C++. Ранее «Рецепты» были только на чистом C и приходилось делать свои классовые обертки. Во-вторых, можно взять одну из АД-библиотек из списка Autodiff.org. По моему опыту использования CPPAD быстрее FADBAD и Trilinos::Sacado примерно на 30%, самая же быстрая, судя по описанию, библиотека это новая ADEPT. В-третьих, для матрично-векторных операций можно взять проверенный временем uBlas , либо новые и быстрые конкуренты Armadillo и blaze-lib — это если для решения СЛАУ использовать отдельные библиотеки (например, SuiteSparce или Pardiso для прямых и ITL для итерационных методов). Привлекательным является также использование интегрированных библиотек линейной алгебры и решателей СЛАУ Eigen, MTL, PETSc (имеются также Ньютоновские решатели на C). Вся триада из заголовка полностью реализована, насколько мне известно, только в Trilinos.
О применении численного дифференцирования
Производные можно вычислять аналитически и численно. К аналитическим методам относятся ручное дифференцирование, символьное (Maple, Wolfram и т.п.) и непосредственно автоматическое дифференцирование, выраженное в средствах выбранного языка программирования.
Современный тренд на использование АД оправдан одной простой причиной — с помощью этой техники устраняется избыточность кода и его дублирование. Другой аргумент состоит в том, что, например, при решении нелинейных дифференциальных уравнений (систем) сеточными методами способ вычисления F(X) сам по себе является нетривиальной задачей. В реальных задачах невязка F(X) представлена суперпозицией вызовов функций с разных слоев программы и ручное дифференцирование теряет свою наглядность. Следует также отметить, что при моделировании нестационарных процессов F(X) меняется на каждом шаге по времени, также может меняться и сам вектор неизвестных X. Использование АД позволяет нам сконцентрироваться непосредственно на формировании F(X), однако не снимает вопрос о верификации получаемой матрицы Якобиана dF(X)/dX. Дело в том, что при вычислении невязок мы можем забыть изменить тип промежуточной переменной со стандартного double на тип АД (а многие библиотеки имеют неявное приведение типа АД к double), тем самым некоторые производные будут вычислены некорректно. Проблема в этом случае состоит в том, что даже при наличии ошибок в формулах для производных метод Ньютона может сходиться, хоть и за возросшее число итераций. Это может быть незаметно при одних начальных данных и приводить к расходимости процесса при других.
Таким образом, какой бы аналитической способ дифференцирования df/dx не был выбран, его крайне желательно дополнить сравнением с численным дифференцированием (f(x + h) — f(x)) / h, иначе всегда будут оставаться сомнения в правильности кода. Естественно, численные производные никогда не совпадут с точностью с правильными аналитическими, тем не менее можно порекомендовать следующий прием юнит-тестирования. Можно выгрузить в текстовые файлы вектора X, F(X) и матрицу dF(X)/dX и выложить на SVN. Затем получить dF(X)/dX численно и сохранить файл поверх существующего, после чего визуально сравнивать файлы между собой. Здесь Вы всегда увидите насколько поменялись значения и сможете локализовать координаты элементов матрицы с большими отклонениями (не в долях) — в этом месте находится ошибка аналитического дифференцирования.
Реализация АД
В Embarcadero Delphi до версии XE5 отсутствует возможность перегрузки арифметических операций для классов, но есть такая возможность для структур record (ссылка):
TAutoDiff = packed record public class operator Equal(a, b: TAutoDiff): Boolean; class operator Negative(v: TAutoDiff): TAutoDiff; class operator Add(a, b: TAutoDiff): TAutoDiff; //и далее по списку end; Обычно в АД на C++ размерность вектора производных является переменной величиной и инициализируется в конструкторе. В Delphi нельзя (есть попытки обойти) перегрузить оператор присваивания для структур и в связи с побитовым копированием вектор производных приходится делать фиксированной длины. Соответствующая константа TAutoDiff.n_all должна изначально задаваться вручную.
Пример 1
procedure TestAutoDiff_1; var i: integer; x, dy: double; x_ad, y_ad: TAutoDiff; begin x := 4; //обязательное зануление всех производных при инициализации x_ad.Init; //сделаем x_ad скалярной константой x_ad := x; assert(x_ad.val = x); for i := 0 to TAutoDiff.n_all - 1 do assert(x_ad.dx[i] = 0); //сравниваем с ручным дифференцированием x_ad.dx[0] := 1; y_ad := sqr(x_ad) + 1 / x_ad; dy := 2 * x - 1 / sqr(x);//вместо x можно использовать x_ad assert(y_ad.dx[0] = dy); end; На данный момент в библиотеке перегружены почти все бинарные и унарные операторы, за исключением тригонометрических функций и булевых операций сдвига. Недостающие функции легко доработать самостоятельно.
Реализация АД для разряженных матриц
С одной стороны фиксированное значение n_all это существенное ограничение, ведь размерность вектора может поступать извне. С другой стороны для некоторых задач его можно ослабить. Дело в том, что если говорить о численных методах решения уравнений механики сплошных сред, то возникающие в них матрицы имеют разреженную структуру — классический пример это «схема крест» для оператора Лапласа, когда в уравнении для узла с координатами (i, j) (ограничимся 2D) задействованы только 5 узлов: (i, j), (i-1, j), (i+1, j), (i, j-1), (i, j+1). Обобщая идею мы должны заложить следующее для данной конкретной задачи:
const n_juncs = 5;//задаем число соседних узлов const n_junc_vars = 1;//задаем число неизвестных в узле const n_all = n_juncs * n_junc_vars; Пусть общее число узлов в решаемой задаче N. Тогда в матрице Якобиана N_all = N * n_junc_vars столбцов, из них ненулевых только n_all. Если завести теперь внутри структуры TAutoDiff целочисленный вектор n_juncs, каждый элемент которого определяет глобальный индекс узла от 0 до N-1, то функцию dx с локальным индексом из диапазона [0, n_all-1] можно дополнить функцией dx_global с уже глобальным индексом из диапазона [0, N_all-1]. Пример 2 иллюстрирует детали использования такого интерфейса, плюсы которого будут наиболее полно видны при реализации метода Ньютона ниже.
Пример 2
procedure TestAutoDiff_2; const N = 100;//размерность расчетной сетки N x N const i = 50; const j = 50; //отображение двумерного индекса узла (i, j) //в построчное хранение в векторе (нумерация с 0) function Splice2d(i, j: integer): integer; begin result := ((i - 1) * N + j - 1); end; var k: integer; n_junc_vars: integer; x: TAutoDiff; juncs_offset: TAutoDiffJuncVector;//вектор begin //n_juncs ячеек с разными смещениями juncs_offset[0] := Splice2d(i - 1, j); juncs_offset[1] := Splice2d(i, j - 1); juncs_offset[2] := Splice2d(i, j); juncs_offset[3] := Splice2d(i, j + 1); juncs_offset[4] := Splice2d(i + 1, j); n_junc_vars := TAutoDiff.n_junc_vars; //задаем, что в векторе dx: //первые n_junc_vars неизвестных относятся к узлу juncs_offset[0] //вторые n_junc_vars неизвестных относятся к узлу juncs_offset[1] //третьи n_junc_vars неизвестных относятся к узлу juncs_offset[2] //четвертые n_junc_vars неизвестных относятся к узлу juncs_offset[3] //последние n_junc_vars неизвестных относятся к узлу juncs_offset[4] x.Init(juncs_offset); //если в dx_global передать индекс вне заданных juncs_offset, //то если он стоит справа от знака равно, то вернет 0, //иначе - сгенерирует исключение for k := 0 to n_junc_vars - 1 do begin x.dx_global[Splice2d(i - 1, j) * n_junc_vars + k] := 1; assert(x.dx[0 * n_junc_vars + k] = 1); x.dx_global[Splice2d(i, j - 1) * n_junc_vars + k] := 1; assert(x.dx[1 * n_junc_vars + k] = 1); x.dx_global[Splice2d(i, j) * n_junc_vars + k] := 1; assert(x.dx[2 * n_junc_vars + k] = 1); x.dx_global[Splice2d(i, j + 1) * n_junc_vars + k] := 1; assert(x.dx[3 * n_junc_vars + k] = 1); x.dx_global[Splice2d(i + 1, j) * n_junc_vars + k] := 1; assert(x.dx[4 * n_junc_vars + k] = 1); end; end; В следующей части планируется рассмотрение класса методов Ньютоновского типа, а также вопроса выбора разряженного решателя СЛАУ.
Пока же предлагаю читателям:
- попробовать написать АД на C++11 с использованием семантики перемещений: 1) это должно работать очень быстро; 2) можно будет обойтись перегрузкой операторов без expression templates; 3) это будет впервые.
- идею для курсовой по реализации АД на Roslyn CTP: можно работать сразу с синтаксическим деревом, которое содержит всю информацию об арифметических выражениях в F(X).
unit AutoDiff; interface const SMALL = 1e-12; const n_juncs = 5; type TAutoDiffJuncVector = array[0..n_juncs - 1] of integer; TAutoDiff = packed record const n_junc_vars = 10; const n_all = n_juncs * n_junc_vars; private juncs_offset: TAutoDiffJuncVector; //<0 если вне диапазона function IsIgnoreJuncOffset: boolean; function IndexOf_dx(glob_indx: integer): integer; function Get_dx_global(glob_indx: integer): double; procedure Set_dx_global(glob_indx: integer; val: double); procedure PrepareBinaryOp(a, b: TAutoDiff); procedure PrepareUnaryOp(v: TAutoDiff); public val: double; dx: array[0..n_all - 1] of double; procedure Init; overload; procedure Init(juncs_offset: TAutoDiffJuncVector); overload; procedure Independent(ir: integer); overload; procedure Independent(juncs_offset: TAutoDiffJuncVector; ir: integer; orient_from_beg: boolean = true); overload; property dx_global[glob_indx: integer]: double read Get_dx_global write Set_dx_global; procedure SetVal(v: TAutoDiff); procedure NoJac(flg: boolean); class operator Implicit(v: double): TAutoDiff; overload; class operator Implicit(v: TAutoDiff): double; overload; class operator Equal(a, b: TAutoDiff): Boolean; class operator Equal(a: double; b: TAutoDiff): Boolean; overload; class operator Equal(a: TAutoDiff; b: double): Boolean; overload; class operator NotEqual(a, b: TAutoDiff): Boolean; class operator NotEqual(a: double; b: TAutoDiff): Boolean; overload; class operator NotEqual(a: TAutoDiff; b: double): Boolean; overload; class operator LessThan(a, b: TAutoDiff): Boolean; class operator LessThan(a: double; b: TAutoDiff): Boolean; overload; class operator LessThan(a: TAutoDiff; b: double): Boolean; overload; class operator LessThanOrEqual(a, b: TAutoDiff): Boolean; class operator LessThanOrEqual(a: double; b: TAutoDiff): Boolean; overload; class operator LessThanOrEqual(a: TAutoDiff; b: double): Boolean; overload; class operator GreaterThan(a, b: TAutoDiff): Boolean; class operator GreaterThan(a: double; b: TAutoDiff): Boolean; overload; class operator GreaterThan(a: TAutoDiff; b: double): Boolean; overload; class operator GreaterThanOrEqual(a, b: TAutoDiff): Boolean; class operator GreaterThanOrEqual(a: double; b: TAutoDiff): Boolean; overload; class operator GreaterThanOrEqual(a: TAutoDiff; b: double): Boolean; overload; class operator Negative(v: TAutoDiff): TAutoDiff; class operator Positive(v: TAutoDiff): TAutoDiff; class operator Add(a, b: TAutoDiff): TAutoDiff; class operator Add(a: double; b: TAutoDiff): TAutoDiff; overload; class operator Add(a: TAutoDiff; b: double): TAutoDiff; overload; class operator Subtract(a, b: TAutoDiff): TAutoDiff; class operator Subtract(a: double; b: TAutoDiff): TAutoDiff; overload; class operator Subtract(a: TAutoDiff; b: double): TAutoDiff; overload; class operator Multiply(a, b: TAutoDiff): TAutoDiff; class operator Multiply(a: double; b: TAutoDiff): TAutoDiff; overload; class operator Multiply(a: TAutoDiff; b: double): TAutoDiff; overload; class operator Divide(a, b: TAutoDiff): TAutoDiff; class operator Divide(a: double; b: TAutoDiff): TAutoDiff; overload; class operator Divide(a: TAutoDiff; b: double): TAutoDiff; overload; end; function sqr(v: TAutoDiff): TAutoDiff; overload; function sqrt(v: TAutoDiff): TAutoDiff; overload; function exp(v: TAutoDiff): TAutoDiff; overload; function ln(v: TAutoDiff): TAutoDiff; overload; //v(x) ^ n(x) function power(a: TAutoDiff; n: double): TAutoDiff; overload; function power(a: double; n: TAutoDiff): TAutoDiff; overload; function power(a: TAutoDiff; n: TAutoDiff): TAutoDiff; overload; function abs(v: TAutoDiff): TAutoDiff; overload; function min(a, b: TAutoDiff): TAutoDiff; function max(a, b: TAutoDiff): TAutoDiff; // function IfThen(flg: Boolean; const on_true: TAutoDiff; const on_false: TAutoDiff): TAutoDiff; function clamp(val, min, max: TAutoDiff): TAutoDiff; //todo: log_a function abs(v: double): double; overload; function sqrt(v: double): double; overload; function sqr(v: double): double; overload; function exp(v: double): double; overload; function ln(v: double): double; overload; procedure swap(var x, y: TAutoDiff); overload; procedure swap(var x, y: double); overload; implementation uses Math, SysUtils; //============================================================================== procedure TAutoDiff.PrepareBinaryOp(a, b: TAutoDiff); var i: integer; begin if a.juncs_offset[0] >= 0 then begin if (b.juncs_offset[0] >= 0) then begin for i := 0 to n_juncs - 1 do begin if a.juncs_offset[i] <> b.juncs_offset[i] then raise Exception.Create('PrepareBinaryOp: must be a.juncs_offset[i] = b.juncs_offset[i]'); end; end; juncs_offset := a.juncs_offset; end else begin juncs_offset := b.juncs_offset; end; end; procedure TAutoDiff.PrepareUnaryOp(v: TAutoDiff); var i: integer; begin juncs_offset := v.juncs_offset;//копирование статического массива (не по ссылке) end; //============================================================================== procedure TAutoDiff.Init; var i: integer; begin for i := 0 to n_juncs - 1 do begin juncs_offset[i] := 0; end; Init(juncs_offset); end; procedure TAutoDiff.Init(juncs_offset: TAutoDiffJuncVector); var i: integer; begin for i := 0 to n_all - 1 do begin dx[i] := 0; end; for i := 0 to n_juncs - 1 do begin self.juncs_offset[i] := juncs_offset[i]; end; end; procedure TAutoDiff.Independent(ir: integer); begin Independent(juncs_offset, ir); end; procedure TAutoDiff.Independent(juncs_offset: TAutoDiffJuncVector; ir: integer; orient_from_beg: boolean); var i, loc_i, junc_i: integer; begin Init(juncs_offset); loc_i := IndexOf_dx(ir); if loc_i >= 0 then begin dx[loc_i] := 1; end else assert(false); end; function TAutoDiff.IsIgnoreJuncOffset: boolean; var i, beg: integer; begin result := true; for i := 1 to n_juncs - 1 do begin if juncs_offset[i] <> juncs_offset[0] then begin result := false; break; end; end; end; function TAutoDiff.IndexOf_dx(glob_indx: integer): integer; var i, offset: integer; begin if IsIgnoreJuncOffset then begin offset := glob_indx - juncs_offset[0] * n_junc_vars; if (0 <= offset) and (offset < n_junc_vars) then result := offset else result := -1; end else begin for i := 0 to n_juncs - 1 do begin offset := glob_indx - juncs_offset[i] * n_junc_vars; if (0 <= offset) and (offset < n_junc_vars) then begin assert(n_junc_vars <= n_junc_vars); if (offset < n_junc_vars) then begin result := i * n_junc_vars + offset; exit; end; end; end; result := -1; end; end; function TAutoDiff.Get_dx_global(glob_indx: integer): double; var loc_i: integer; begin loc_i := IndexOf_dx(glob_indx); if loc_i >= 0 then result := dx[loc_i] else result := 0; end; procedure TAutoDiff.Set_dx_global(glob_indx: integer; val: double); var loc_i: integer; begin loc_i := IndexOf_dx(glob_indx); if loc_i >= 0 then dx[loc_i] := val else assert(false); end; procedure TAutoDiff.SetVal(v: TAutoDiff); begin val := v.val; end; class operator TAutoDiff.Implicit(v: double): TAutoDiff; begin result.val := v; result.NoJac(true); end; procedure TAutoDiff.NoJac(flg: boolean); const NO_JAC_MARK = -1; var i: integer; begin Init; if flg then begin for i := 0 to n_juncs - 1 do begin juncs_offset[i] := NO_JAC_MARK; end; end; end; class operator TAutoDiff.Implicit(v: TAutoDiff): double; begin result := v.val; end; class operator TAutoDiff.Equal(a, b: TAutoDiff): Boolean; begin result := (a.val = b.val); end; class operator TAutoDiff.Equal(a: double; b: TAutoDiff): Boolean; begin result := (a = b.val); end; class operator TAutoDiff.Equal(a: TAutoDiff; b: double): Boolean; begin result := (a.val = b); end; //------------------------------------------------------------------------------ class operator TAutoDiff.NotEqual(a, b: TAutoDiff): Boolean; begin result := (a.val <> b.val); end; class operator TAutoDiff.NotEqual(a: double; b: TAutoDiff): Boolean; begin result := (a <> b.val); end; class operator TAutoDiff.NotEqual(a: TAutoDiff; b: double): Boolean; begin result := (a.val <> b); end; //------------------------------------------------------------------------------ class operator TAutoDiff.LessThan(a, b: TAutoDiff): Boolean; begin result := (a.val < b.val); end; class operator TAutoDiff.LessThan(a: double; b: TAutoDiff): Boolean; begin result := (a < b.val); end; class operator TAutoDiff.LessThan(a: TAutoDiff; b: double): Boolean; begin result := (a.val < b); end; //------------------------------------------------------------------------------ class operator TAutoDiff.LessThanOrEqual(a, b: TAutoDiff): Boolean; begin result := (a.val <= b.val); end; class operator TAutoDiff.LessThanOrEqual(a: double; b: TAutoDiff): Boolean; begin result := (a <= b.val); end; class operator TAutoDiff.LessThanOrEqual(a: TAutoDiff; b: double): Boolean; begin result := (a.val <= b); end; //------------------------------------------------------------------------------ class operator TAutoDiff.GreaterThan(a, b: TAutoDiff): Boolean; begin result := (a.val > b.val); end; class operator TAutoDiff.GreaterThan(a: double; b: TAutoDiff): Boolean; begin result := (a > b.val); end; class operator TAutoDiff.GreaterThan(a: TAutoDiff; b: double): Boolean; begin result := (a.val > b); end; //------------------------------------------------------------------------------ class operator TAutoDiff.GreaterThanOrEqual(a, b: TAutoDiff): Boolean; begin result := (a.val >= b.val); end; class operator TAutoDiff.GreaterThanOrEqual(a: double; b: TAutoDiff): Boolean; begin result := (a >= b.val); end; class operator TAutoDiff.GreaterThanOrEqual(a: TAutoDiff; b: double): Boolean; begin result := (a.val >= b); end; //------------------------------------------------------------------------------ class operator TAutoDiff.Negative(v: TAutoDiff): TAutoDiff; begin result := - 1 * v; end; class operator TAutoDiff.Positive(v: TAutoDiff): TAutoDiff; begin result := v; end; //------------------------------------------------------------------------------ class operator TAutoDiff.Add(a, b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := a.val + b.val; result.PrepareBinaryOp(a, b); for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i] + b.dx[i]; end; end; class operator TAutoDiff.Add(a: double; b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := a + b.val; result.PrepareUnaryOp(b); for i := 0 to result.n_all - 1 do begin result.dx[i] := b.dx[i]; end; end; class operator TAutoDiff.Add(a: TAutoDiff; b: double): TAutoDiff; var i: integer; begin result.val := a.val + b; result.PrepareUnaryOp(a); for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i]; end; end; //------------------------------------------------------------------------------ class operator TAutoDiff.Subtract(a, b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := a.val - b.val; result.PrepareBinaryOp(a, b); for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i] - b.dx[i]; end; end; class operator TAutoDiff.Subtract(a: double; b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := a - b.val; result.PrepareUnaryOp(b); for i := 0 to result.n_all - 1 do begin result.dx[i] := - b.dx[i]; end; end; class operator TAutoDiff.Subtract(a: TAutoDiff; b: double): TAutoDiff; var i: integer; begin result.val := a.val - b; result.PrepareUnaryOp(a); for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i]; end; end; //------------------------------------------------------------------------------ class operator TAutoDiff.Multiply(a, b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := a.val * b.val; result.PrepareBinaryOp(a, b); for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i] * b.val + a.val * b.dx[i]; end; end; class operator TAutoDiff.Multiply(a: double; b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := a * b.val; result.PrepareUnaryOp(b); for i := 0 to result.n_all - 1 do begin result.dx[i] := a * b.dx[i]; end; end; class operator TAutoDiff.Multiply(a: TAutoDiff; b: double): TAutoDiff; var i: integer; begin result.val := a.val * b; result.PrepareUnaryOp(a); for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i] * b; end; end; //------------------------------------------------------------------------------ class operator TAutoDiff.Divide(a, b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := a.val / b.val; result.PrepareBinaryOp(a, b); for i := 0 to result.n_all - 1 do begin result.dx[i] := (a.dx[i] * b.val - a.val * b.dx[i]) / System.sqr(b.val); end; end; class operator TAutoDiff.Divide(a: double; b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := a / b.val; result.PrepareUnaryOp(b); for i := 0 to result.n_all - 1 do begin result.dx[i] := - a * b.dx[i] / System.sqr(b.val); end; end; class operator TAutoDiff.Divide(a: TAutoDiff; b: double): TAutoDiff; var i: integer; begin result.val := a.val / b; result.PrepareUnaryOp(a); for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i] / b; end; end; //============================================================================== function sqr(v: TAutoDiff): TAutoDiff; var i: integer; d: double; begin result.val := System.sqr(v.val); result.PrepareUnaryOp(v); d := 2 * v.val; for i := 0 to v.n_all - 1 do begin result.dx[i] := d * v.dx[i]; end; end; function sqrt(v: TAutoDiff): TAutoDiff; var i: integer; d: double; begin result.val := System.sqrt(v.val); result.PrepareUnaryOp(v); if abs(result.val) < SMALL then d := 0 else d := 0.5 / result.val; for i := 0 to v.n_all - 1 do begin result.dx[i] := d * v.dx[i]; end; end; function exp(v: TAutoDiff): TAutoDiff; var i: integer; begin result.val := System.exp(v.val); result.PrepareUnaryOp(v); for i := 0 to v.n_all - 1 do begin result.dx[i] := result.val * v.dx[i]; end; end; function ln(v: TAutoDiff): TAutoDiff; var i: integer; begin result.val := System.ln(v.val); result.PrepareUnaryOp(v); for i := 0 to v.n_all - 1 do begin result.dx[i] := v.dx[i] / v.val; end; end; function power(a: TAutoDiff; n: double): TAutoDiff; var i: integer; d: double; begin result.val := Math.power(a.val, n); result.PrepareUnaryOp(a); d := n * Math.power(a.val, n - 1); for i := 0 to result.n_all - 1 do begin result.dx[i] := d * a.dx[i]; end; end; function power(a: double; n: TAutoDiff): TAutoDiff; var i: integer; d: double; begin result.val := Math.power(a, n.val); result.PrepareUnaryOp(n); d := ln(n) * result.val; for i := 0 to n.n_all - 1 do begin result.dx[i] := d * n.dx[i]; end; end; function power(a: TAutoDiff; n: TAutoDiff): TAutoDiff; var i: integer; d: double; begin result.val := Math.power(a.val, n.val); result.PrepareUnaryOp(n); d := Math.power(a.val, n.val - 1); for i := 0 to n.n_all - 1 do begin result.dx[i] := d * (n.val * a.dx[i] + a.val * ln(a.val) * n.dx[i]); end; end; function abs(v: TAutoDiff): TAutoDiff; var i: integer; begin result.val := System.abs(v.val); result.PrepareUnaryOp(v); for i := 0 to v.n_all - 1 do begin result.dx[i] := Math.sign(v.val) * v.dx[i]; end; end; function min(a, b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := Math.min(a.val, b.val); result.PrepareBinaryOp(a, b); if a.val < b.val then begin for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i]; end; end else begin for i := 0 to result.n_all - 1 do begin result.dx[i] := b.dx[i]; end; end; end; function max(a, b: TAutoDiff): TAutoDiff; var i: integer; begin result.val := Math.max(a.val, b.val); result.PrepareBinaryOp(a, b); if a.val > b.val then begin for i := 0 to result.n_all - 1 do begin result.dx[i] := a.dx[i]; end; end else begin for i := 0 to result.n_all - 1 do begin result.dx[i] := b.dx[i]; end; end; end; function IfThen(flg: Boolean; const on_true: TAutoDiff; const on_false: TAutoDiff): TAutoDiff; begin if flg then result := on_true else result := on_false; end; function clamp(val, min, max: TAutoDiff): TAutoDiff; begin Result := IfThen(val < min, min, IfThen(max < val, max, val)); end; procedure Swap(var x, y: TAutoDiff); var tmp: TAutoDiff; begin tmp := x; x := y; y := tmp; end; procedure Swap(var x, y: double); var tmp: double; begin tmp := x; x := y; y := tmp; end; //============================================================================== function abs(v: double): double; begin result := System.Abs(v); end; function sqrt(v: double): double; begin result := System.sqrt(v); end; function sqr(v: double): double; begin result := System.sqr(v); end; function exp(v: double): double; begin result := system.Exp(v); end; function ln(v: double): double; begin result := system.Ln(v); end; end.
ссылка на оригинал статьи http://habrahabr.ru/post/247379/
Добавить комментарий