Delphiの選択はレガシーコードによって決定されますが、C ++では次のように問題を解決できることを最初に注意しましょう。 まず、このタイプのニュートン法（基本的なニュートン法、ブロイデン法）には、 C ++の数値レシピがあります 。 以前は、レシピは純粋なCのみであり、独自のクラスラッパーを作成する必要がありました。 次に、 Autodiff.orgリストからADライブラリの1つを取得できます。 私の経験では、CPPADの使用はFADBADおよびTrilinos :: Sacadoよりも約30％高速ですが、説明から判断すると、最も速いのは新しいADEPTライブラリです。 第三に、行列ベクトル演算では、 SLAEを解決するために個別のライブラリを使用する場合（たとえば、直接法にはSuiteSparceまたはPardiso 、反復法にはITL ）、 uBlas 、時間テスト済み 、または新しい高速のArmadilloおよびblaze-libを使用できます。 統合線形代数ライブラリとSLAE Eigen、 MTL 、 PETScのソルバーを使用することも魅力的です （ Cにはニュートンのソルバーもあります）。 私の知る限り、ヘッダーからのトライアド全体は、 Trilinosでのみ完全に実装されています 。
数値微分の適用について
デリバティブは、分析的および数値的に計算できます。 分析方法には、選択したプログラミング言語の手段で表現された、手動による差別化、シンボリック（Maple、Wolframなど）および直接的な自動差別化が含まれます。
血圧の使用に関する現在の傾向は、1つの単純な理由によって正当化されます-この手法を使用すると、コードの冗長性と重複が排除されます。 もう1つの議論は、たとえば、グリッド法によって非線形微分方程式（システム）を解くとき、F（X）を計算する方法自体は重要なタスクであるということです。 実際の問題では、不一致F（X）は異なるプログラムレイヤーからの関数呼び出しの重ね合わせによって表され、手動による区別はその可視性を失います。 また、非定常プロセスをモデル化する場合、F（X）が各タイムステップで変化するため、未知のXのベクトル自体も変化する可能性があることに注意してください。ADを使用すると、F（X）ヤコビアンdF（X）/ dX。 事実、残差を計算するとき、中間変数の型を標準のdoubleからHELLの型に変更することを忘れることがあります（そして多くのライブラリーはHELL型の暗黙的なキャストをdoubleに持っています）。 この場合の問題は、微分の式に誤差があったとしても、反復回数が増えてもニュートン法が収束する可能性があることです。 これは、一部の初期データでは知覚できない場合があり、他のプロセスとのプロセスの相違につながる可能性があります。
したがって、どのような分析の微分df / dxが選択されたとしても、数値微分（f（x + h）-f（x））/ hとの比較で補完することが非常に望ましいです。そうでなければ、コードの正確性について常に疑問が生じます。 当然、数値微分は正しい分析の精度と一致することはありませんが、それでも次の単体テスト手法を推奨できます。 X、F（X）ベクトルおよびdF（X）/ dXマトリックスをテキストファイルにアップロードし、SVNに配置できます。 次に、dF（X）/ dXを数値で取得し、既存のファイルの上にファイルを保存してから、ファイルを相互に視覚的に比較します。 ここでは、値がどのように変化したかを常に確認でき、大きな偏差（分数ではない）でマトリックス要素の座標をローカライズできます-これが分析微分の誤差です。
BPの実装
XE5より前のEmbarcadero Delphiでは、クラスの算術演算をオーバーロードする方法はありませんが、レコード構造（リンク）にはこのような可能性があります。
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 ++ ADでは、導関数ベクトルの次元は変数であり、コンストラクターで初期化されます。 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;
現時点では、三角関数とブールシフト演算を除き、ほとんどすべての二項演算子と単項演算子がライブラリにオーバーロードされています。 欠落している関数は、個別に簡単に変更できます。
スパース行列のBP実装
一方では、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に定義する場合、範囲[0、n_all-1]のローカルインデックスを持つdx関数は、範囲[ 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;
次のパートでは、ニュートン型のメソッドのクラスと、スパースSLAEソルバーを選択する問題について検討する予定です。 それまでの間、私は読者に以下を提供します。
- 動作のセマンティクスを使用して、C ++ 11でADを記述してみてください。1）非常に高速に動作するはずです。 2）式テンプレートなしで演算子のオーバーロードを行うことができます。 3）これは初めてです。
- Roslyn CTPでのADの実装に関するコースのアイデア：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.