On automatic differentiation, Newton's method and solution of SLAE on Delphi. Part 1

About automatic differentiation (AD) on Habré has already been written here and here . This article proposes the implementation of AD for Delphi (tested in Embarcadero XE2, XE6) along with convenient classes of Newton methods for solving nonlinear equations f (x) = 0 and systems F (X) = 0. Any references to ready-made similar libraries are welcome, I did not find a similar one, not counting the excellent sparse matrix SLAE solver (see under the cat).

Let's note at the very beginning that the choice of Delphi is due to the legacy code; nevertheless, in C ++, the problem can be solved as follows. First, for Newton's methods (Newton's base method, Broyden method) there are Numerical Recipes in C ++ . Previously, "Recipes" were only on pure C and had to make their own class wrappers. Secondly, you can take one of the AD libraries from the list of Autodiff.org . In my experience using CPPAD faster than FADBAD and Trilinos :: Sacado by about 30%, the fastest, judging by the description, the library is the new ADEPT . Thirdly, it is possible to take the time-tested uBlas for matrix-vector operations, or Armadillo and blaze-lib's new and fast competitors - if separate libraries are used to solve the SLAE (for example, SuiteSparce or Pardiso for direct and ITL for iterative methods). The use of integrated libraries of linear algebra and solvers of the SLAE Eigen , MTL , PETSc is also attractive (there are also Newtonian solvers in C ). The whole triad from the header is fully implemented, as far as I know, only in Trilinos .

On the application of numerical differentiation

Derivatives can be calculated analytically and numerically. Analytical methods include manual differentiation, symbolic (Maple, Wolfram, etc.) and automatic differentiation itself, expressed in the means of the selected programming language.

The current trend in the use of blood pressure is justified by one simple reason - using this technique eliminates the redundancy of the code and its duplication. Another argument is that, for example, when solving nonlinear differential equations (systems) by grid methods, the method of calculating F (X) is in itself a non-trivial task . In real problems, the residual F (X) is represented by a superposition of function calls from different layers of the program, and manual differentiation loses its visibility. It should also be noted that when modeling non-stationary processes, F (X) changes at each time step, the vector of unknowns X can also change. The use of AD allows us to concentrate directly on the formation of F (X), but does not remove the question of the verification of the resulting matrix Jacobian dF (X) / dX. The fact is that when calculating discrepancies, we can forget to change the type of intermediate variable from standard double to type of blood pressure (and many libraries have an implicit type conversion to double), thus some derivatives will be incorrectly calculated. The problem in this case is that even if there are errors in the formulas for the derivatives, Newton's method may converge, albeit for an increased number of iterations. This may be imperceptible with some initial data and lead to a divergence of the process with others.
')
Thus, whatever analytical method of df / dx differentiation is chosen, it is highly desirable to supplement it with a comparison with numerical differentiation (f (x + h) - f (x)) / h, otherwise doubts about the correctness of the code will always remain. Naturally, numerical derivatives will never coincide with accuracy with the correct analytical ones; nevertheless, we can recommend the following unit-testing technique. You can upload the vector files X, F (X) and the matrix dF (X) / dX into text files and upload them to SVN. Then get dF (X) / dX numerically and save the file over the existing one, then visually compare the files with each other. Here you will always see how much the values have changed and you can localize the coordinates of the elements of the matrix with large deviations (not in fractions) - in this place is the error of analytical differentiation.

Implementation of blood pressure

Up to XE5 version of Embarcadero Delphi there is no possibility of overloading arithmetic operations for classes, but there is such an opportunity for record structures (link) :

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;

Usually in AD on C ++, the dimension of the vector of derivatives is variable and initialized in the constructor. In Delphi it is impossible ( there are attempts to bypass ) to overload the assignment operator for structures and, in connection with bitwise copying, the vector of derivatives has to be made of a fixed length. The corresponding TAutoDiff.n_all constant must be initially set manually.

Example 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;

At the moment, almost all binary and unary operators are overloaded in the library, with the exception of trigonometric functions and Boolean shift operations. Missing features are easy to customize.

Implementation of blood pressure for sparse matrices

On the one hand, the fixed value n_all is a significant limitation, because the dimension of the vector can come from outside. On the other hand, for some tasks it can be weakened. The fact is that if we talk about numerical methods for solving equations of continuum mechanics, then the matrices arising in them have sparse structure - a classic example is the “cross scheme” for the Laplace operator, when in the equation for a node with coordinates (i, j) 2D) only 5 nodes are involved: (i, j), (i-1, j), (i + 1, j), (i, j-1), (i, j + 1). Summarizing the idea, we should lay down the following for this particular task:

 const n_juncs = 5;//    const n_junc_vars = 1;//     const n_all = n_juncs * n_junc_vars;

Let the total number of nodes in the problem being solved N. Then, in the Jacobian matrix, N_all = N * n_junc_vars columns, of which non-zero are only n_all. If we now have an integer vector n_juncs inside the structure of TAutoDiff, each element of which defines the global node index from 0 to N-1, then the dx function with a local index from the range [0, n_all-1] can be added with the function dx_global with a global index from the range [ 0, N_all-1]. Example 2 illustrates the details of using such an interface, the advantages of which will be most clearly seen when implementing the Newton method below.

Example 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;

In the next part we plan to consider the class of Newton type methods, as well as the question of choosing a sparse solver of SLAE. In the meantime, I offer readers:

try to write AD on C ++ 11 using displacement semantics: 1) it should work very quickly; 2) it will be possible to do with operator overloading without expression templates; 3) it will be the first time.
idea for coursework on the implementation of blood pressure on Roslyn CTP: you can work immediately with the syntax tree, which contains all the information about arithmetic expressions in F (X).

AutoDiff.pas library file

 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.

Source: https://habr.com/ru/post/247379/

All Articles