Об автоматическом дифференцировании, методе Ньютона и решении СЛАУ на Delphi. Часть 1 / forpes.ru

Главная
Об автоматическом дифференцировании, методе Ньютона и решении СЛАУ на Delphi. Часть 1

Об автоматическом дифференцировании, методе Ньютона и решении СЛАУ на Delphi. Часть 1 +6

09.05.2015 16:25

shtr 9 6832 Источник

Об автоматическом дифференцировании (АД) на Хабре уже писалось здесь и здесь. В данной статье предлагается реализация АД для 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).

Файл библиотеки AutoDiff.pas

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.

Комментарии (9)

docker1
09.05.2015 21:07
#8411399
Матрицы разреженные, а не разряженные все-таки
1. shtr Автор
  09.05.2015 21:28
  #8411413
  Спасибо)

lany
10.05.2015 14:01
#8411645
+2
А почему ненормальное программирование? Или в наши дни любая программа на Delphi считается ненормальным программированием? :-)
1. shtr Автор
  10.05.2015 18:50
  #8411751
  Решил поддержать авторов первых двух ссылок
1. MaxFactor
  11.05.2015 22:40
  #8412623
  Я бы сказал, да. Тут недавно сделал простенькую программу, ничего страшного не делала, а антивирус ругается. Я залез на тоталвирус и загрузил программульку — 16 антивирусов показали, что там вирус, обидно., Как я прочитал на форумах, теперь, чтобы твоя Delphi программа не считалась подозрительной надо сертификат разработчика ставить.
  1. MacIn
    12.05.2015 01:47
    #8412755
    Не встречал такого, хотя наша клиентская база — не единицы и не тысячи. С другой стороны, у меня на машине антивирус не пускает Hello world на С, собранный на VC2010.
  1. proto
    12.05.2015 11:25
    #8413121
    надо просто чтобы манифест был, сам встречался с такой проблемой.
    
    MaxFactor
    12.05.2015 12:15
    #8413241
    Пробовал, согласен получше, но все равно проблема остается. Как я понял, на текущий момент антивирусы всю работу стараются сложить на «сертификацию».

maybe_im_a_leo
12.05.2015 02:45
#8412765
Нужно еще понимать, что не каждую систему нелинейных уравнений можно решить с помощью метода Ньютона. Проблемы начинаются если где-то на траектории матрица Якоби dF(X)/dX вырождается.

Надо либо про регуляризацию думать, либо использовать более хардкорные методы типа продолжения решения по параметру (homotopy parameter continuation)

Об автоматическом дифференцировании, методе Ньютона и решении СЛАУ на Delphi. Часть 1 +6

О применении численного дифференцирования

Реализация АД

Пример 1

Реализация АД для разреженных матриц

Пример 2

Комментарии (9)

docker1

shtr Автор

lany

shtr Автор

MaxFactor

MacIn

proto

MaxFactor

maybe_im_a_leo