const
CONST_1_DIV_SQRT_2 = 0.70710678118654752440084436210485;
FIND_MIN_OK
= 0;
FIND_MIN_INVALID_OPTION = 1;
FIND_MIN_INVALID_FUNC = 2;
FIND_MIN_INVALID_N = 3;
FIND_MIN_INVALID_X0 = 4;
FIND_MIN_INVALID_X = 5;
FIND_MIN_INVALID_EPS = 6;
FIND_MIN_INVALID_DELTA = 7;
FIND_MIN_INVALID_R = 8;
FIND_MIN_MODE_NOT_SUPPORT = 9;
FIND_MIN_OUT_OF_MEMORY = 10;
FIND_MIN_INVALID_ALPHA
= 11;
FIND_MIN_INVALID_BETA = 12;
FIND_MIN_INVALID_GAMMA = 13;
FIND_MIN_MODE_STD = 0;
FIND_MIN_MODE_1 = 1;
FIND_MIN_MODE_2 = 2;
FIND_MIN_USE_EPS = $00000001;
FIND_MIN_USE_R = $00000002;
FIND_MIN_USE_MODE = $00000004;
FIND_MIN_USE_DELTA = $00000008;
FIND_MIN_USE_ALPHA = $00000010;
FIND_MIN_USE_BETA = $00000020;
FIND_MIN_USE_GAMMA = $00000040;
//
Некоторые комбинации стандартных опций:
FIND_MIN_USE_EPS_R = FIND_MIN_USE_EPS or FIND_MIN_USE_R;
FIND_MIN_USE_EPS_R_MODE = FIND_MIN_USE_EPS_R or
FIND_MIN_USE_MODE;
FIND_MIN_USE_EPS_R_DELTA = FIND_MIN_USE_EPS_R or
FIND_MIN_USE_DELTA;
FIND_MIN_USE_EPS_R_MODE_DELTA = FIND_MIN_USE_EPS_R_MODE or
FIND_MIN_USE_DELTA;
FIND_MIN_USE_ALL = FIND_MIN_USE_EPS or FIND_MIN_USE_R or
FIND_MIN_USE_MODE or
FIND_MIN_USE_DELTA or FIND_MIN_USE_ALPHA or
FIND_MIN_USE_BETA or FIND_MIN_USE_GAMMA;
type
PMathFunction = ^TMathFunction;
TMathFunction = function(X: PExtended): Extended;
PNelderOption = ^TNelderOption;
TNelderOption
= record
Size:
Cardinal; // Размер структуры (обязательно)
Flags:
Cardinal; // Флаги (обязательно)
Func:
TMathFunction; // Функция (обязательно)
N:
Integer; // Размерность (обязательно)
X0:
PExtended; // Указатель на начальную точку (обязательно)
X:
PExtended; // Указатель куда записывать результат (обязательно)
Eps:
Extended; // Точность (опция FIND_MIN_USE_EPS)
Delta:
Extended; // Способ проверки (опция FIND_MIN_USE_DELTA)
R:
Extended; // Расстояние между вершинами симплекса (опция FIND_MIN_USE_R)
Mode:
Integer; // Метод решения (опция FIND_MIN_USE_MODE)
Alpha:
Extended; // Коэффициент отражения (опция FIND_MIN_USE_ALPHA)
Beta:
Extended; // Коэффициент сжатия (опция FIND_MIN_USE_BETA)
Gamma:
Extended; // Коэффициент растяжения (опция FIND_MIN_USE_GAMMA)
end;
{**********
Проверка
указателя Option на то, что все его параметры доступны для чтения
**********}
function
CheckNelderOptionPtr(Option: PNelderOption): Integer;
begin
//
Проверка указателя #1 (допустимость указателя)
if
IsBadReadPtr(@Option, 4) then
begin
Result :=
FIND_MIN_INVALID_OPTION;
Exit;
end;
//
Проверка указателя #2 (слишком мало параметров)
if
Option.Size < 24 then
begin
Result := FIND_MIN_INVALID_OPTION;
Exit;
end;
//
Проверка указателя #3 (все данные можно читать?)
if
IsBadReadPtr(@Option, Option.Size) then
begin
Result :=
FIND_MIN_INVALID_OPTION;
Exit;
end;
Result :=
FIND_MIN_OK;
end;
{************
Копирование
данных из одной структуры в другую с попутной проверкой
на
допустимость значений всех параметров.
************}
function
CopyData(const InOption: TNelderOption; var OutOption:
TNelderOption):
Integer;
var
CopyCount:
Cardinal;
begin
Result :=
FIND_MIN_OK;
//
Копируем одну структуру в другую
CopyCount
:= SizeOf(TNelderOption);
if
InOption.Size < CopyCount then
CopyCount
:= InOption.Size;
Move(InOption, OutOption, CopyCount);
//
Устанавливаем размер
OutOption.Size
:= SizeOf(TNelderOption);
//
Проверка Option.Func
if
IsBadCodePtr(@OutOption.Func) then
begin
Result :=
FIND_MIN_INVALID_FUNC;
Exit;
end;
//
Проверка Option.N
if
OutOption.N <= 0 then
begin
Result :=
FIND_MIN_INVALID_N;
Exit;
end;
//
Проверка Option.X0
if
IsBadReadPtr(OutOption.X0, OutOption.N * SizeOf(Extended)) then
begin
Result :=
FIND_MIN_INVALID_X0;
Exit;
end;
//
Проверка Option.X
if
IsBadWritePtr(OutOption.X, OutOption.N * SizeOf(Extended)) then
begin
Result :=
FIND_MIN_INVALID_X;
Exit;
end;
//
Проверка Option.Eps
if
(FIND_MIN_USE_EPS and OutOption.Flags) <> 0 then
begin
if
OutOption.Size < 34 then // Eps не вписывается в размер
begin
Result
:= FIND_MIN_INVALID_OPTION;
Exit;
end
else
if OutOption.Eps <= 0 then
begin
Result
:= FIND_MIN_INVALID_EPS;
Exit;
end;
end
else
begin
OutOption.Eps := 1E-06; // Default value;
end;
//
Проверка OutOption.Delta
if (FIND_MIN_USE_DELTA
and OutOption.Flags) <> 0 then
begin
if
OutOption.Size < 44 then
begin
Result
:= FIND_MIN_INVALID_OPTION;
Exit;
end
else
if (OutOption.Delta < 0.0) or (OutOption.Delta > 1.0) then
begin
Result
:= FIND_MIN_INVALID_DELTA;
Exit;
end;
end
else
begin
OutOption.Delta := 0.5; // Default value
end;
//
Проверка OutOption.R
if
(FIND_MIN_USE_R and OutOption.Flags) <> 0 then
begin
if
OutOption.Size < 54 then
begin
Result
:= FIND_MIN_INVALID_OPTION;
Exit;
end
else
if (OutOption.R <= 0.0) then
begin
Result
:= FIND_MIN_INVALID_R;
Exit;
end;
end
else
begin
OutOption.R := 100.0; // Default value
end;
//
Проверка OutOption.Mode
if
(FIND_MIN_USE_MODE and OutOption.Flags) <> 0 then
begin
if
OutOption.Size < 58 then
begin
Result
:= FIND_MIN_INVALID_OPTION;
Exit;
end
else
if (OutOption.Mode <> FIND_MIN_MODE_STD) then
if
(OutOption.Mode <> FIND_MIN_MODE_1) then
if
(OutOption.Mode <> FIND_MIN_MODE_2) then
begin
Result := FIND_MIN_MODE_NOT_SUPPORT;
Exit;
end;
end
else
begin
OutOption.Mode := FIND_MIN_MODE_STD; // Default value
end;
//
Проверка OutOption.Alpha
if
(FIND_MIN_USE_ALPHA and OutOption.Flags) <> 0 then
begin
if
OutOption.Size < 68 then
begin
Result
:= FIND_MIN_INVALID_OPTION;
Exit;
end
else
if OutOption.Alpha < 0.0 then
begin
Result
:= FIND_MIN_INVALID_ALPHA;
Exit;
end;
end
else
begin
OutOption.Alpha := 1.0; // Default value
end;
//
Проверка OutOption.Beta
if
(FIND_MIN_USE_BETA and OutOption.Flags) <> 0 then
begin
if
OutOption.Size < 78 then
begin
Result
:= FIND_MIN_INVALID_OPTION;
Exit;
end
else
if (OutOption.Beta < 0.0) or (OutOption.Beta > 1.0) then
begin
Result
:= FIND_MIN_INVALID_BETA;
Exit;
end;
end
else
begin
OutOption.Beta := 0.5; // Default value
end;
//
Проверка OutOption.Gamma
if
(FIND_MIN_USE_GAMMA and OutOption.Flags) <> 0 then
begin
if
OutOption.Size < 88 then
begin
Result
:= FIND_MIN_INVALID_OPTION;
Exit;
end
else
if OutOption.Gamma < 1.0 then
begin
Result
:= FIND_MIN_INVALID_GAMMA;
Exit;
end;
end
else
begin
OutOption.Gamma := 1.5; // Default value
end;
end;
type
TNelderTempData = record
D1:
Extended;
D2:
Extended;
FC:
Extended;
FU:
Extended;
XN:
PExtended;
D0:
PExtended;
FX:
PExtended;
C:
PExtended;
U:
PExtended;
V:
PEXtended;
Indexes:
PInteger;
end;
function
InitializeTempData(var TempData: TNelderTempData; N: Integer): Integer;
var
SizeD0:
Integer;
SizeFX:
Integer;
SizeC:
Integer;
SizeU:
Integer;
SizeV:
Integer;
SizeIndexes: Integer;
SizeAll:
Integer;
Ptr: PChar;
begin
//
Вычисляем размеры
SizeD0 := N
* (N + 1) * SizeOf(Extended);
SizeFX :=
(N + 1) * SizeOf(Extended);
SizeC := N
* SizeOf(Extended);
SizeU := N
* SizeOf(Extended);
SizeV := N
* SizeOf(Extended);
SizeIndexes
:= (N + 1) * SizeOf(Integer);
SizeAll :=
SizeD0 + SizeFX + SizeC + SizeU + SizeV + SizeIndexes;
Ptr :=
SysGetMem(SizeAll);
if not
Assigned(Ptr) then
begin
Result :=
FIND_MIN_OUT_OF_MEMORY;
Exit;
end;
TempData.D0
:= PExtended(Ptr);
Ptr := Ptr
+ SizeD0;
TempData.FX
:= PExtended(Ptr);
Ptr := Ptr
+ SizeFX;
TempData.C
:= PExtended(Ptr);
Ptr := Ptr
+ SizeC;
TempData.U
:= PExtended(Ptr);
Ptr := Ptr
+ SizeU;
TempData.V
:= PExtended(Ptr);
Ptr := Ptr
+ SizeV;
TempData.Indexes := PInteger(Ptr);
// Ptr
:= Ptr + SizeIndexes
Result :=
FIND_MIN_OK;
end;
procedure
FinalizeTempData(var TempData: TNelderTempData);
var
Ptr:
Pointer;
begin
Ptr :=
TempData.D0;
TempData.D0
:= nil;
TempData.FX
:= nil;
TempData.C
:= nil;
TempData.U
:= nil;
TempData.V
:= nil;
TempData.Indexes := nil;
SysFreeMem(Ptr);
end;
{*********
Строится
симплекс:
В целях
оптимизации поменяем местами строки и столбцы
**********}
procedure
BuildSimplex(var Temp: TNelderTempData; const Option:
TNelderOption);
var
I, J:
Integer;
PtrD0:
PExtended;
begin
with
Temp, Option do
begin
D1 :=
CONST_1_DIV_SQRT_2 * R * (Sqrt(N + 1) + N - 1) / N;
D2 :=
CONST_1_DIV_SQRT_2 * R * (Sqrt(N + 1) - 1) / N;
FillChar(D0^, N * SizeOf(Extended), 0);
PtrD0 :=
D0;
PChar(PtrD0) := PChar(PtrD0) + N * SizeOf(Extended);
for
I := 0 to N - 1 do
for
J := 0 to N - 1 do
begin
if
I = J then
PtrD0^ := D1
else
PtrD0^ := D2;
PChar(PtrD0) := PChar(PtrD0) + SizeOf(Extended);
end;
end;
end;
{*********
Перемещение
симплекса в точку A
*********}
procedure
MoveSimplex(var Temp: TNelderTempData; const Option:
TNelderOption; A:
PExtended);
var
I, J:
Integer;
PtrA,
PtrD0: PExtended;
begin
with
Temp, Option do
begin
PtrD0 :=
D0;
for
I := 0 to N do
begin
PtrA :=
A;
for
J := 0 to N - 1 do
begin
PtrD0^ := PtrD0^ + PtrA^;
PChar(PtrD0) := PChar(PtrD0) + SizeOf(Extended);
PChar(PtrA) := PChar(PtrA) + SizeOf(Extended);
end;
end;
end;
end;
// Быстрая сортировка значений FX
procedure
QuickSortFX(L, R: Integer; FX: PExtended; Indexes: PInteger);
var
I, J, K:
Integer;
P, T:
Extended;
begin
repeat
I := L;
J := R;
// P
:= FX[(L+R) shr 1] :
P :=
PExtended(PChar(FX) + SizeOf(Extended) * ((L + R) shr 1))^;
repeat
//
while FX[I] < P do Inc(I):
while
PExtended(PChar(FX) + SizeOf(Extended) * I)^ < P do
Inc(I);
//
while FX[J] > P do Dec(J):
while
PExtended(PChar(FX) + SizeOf(Extended) * J)^ > P do
Dec(J);
if
I <= J then
begin
//
Переставляем местами значения FX
T :=
PExtended(PChar(FX) + SizeOf(Extended) * I)^;
PExtended(PChar(FX) + SizeOf(Extended) * I)^ := PExtended(PChar(FX) +
SizeOf(Extended) * J)^;
PExtended(PChar(FX) + SizeOf(Extended) * J)^ := T;
//
Поддерживаем порядок и в индексах
K :=
PInteger(PChar(Indexes) + SizeOf(Integer) * I)^;
PInteger(PChar(Indexes)
+ SizeOf(Integer) * I)^ :=
PInteger(PChar(Indexes) + SizeOf(Integer) * J)^;
PInteger(PChar(Indexes) + SizeOf(Integer) * J)^ := K;
Inc(I);
Dec(J);
end;
until
I > J;
if
L < J then
QuickSortFX(L, J, FX, Indexes);
L := I;
until
I >= R;
end;
procedure
CalcFX(var Temp: TNelderTempData; Option: TNelderOption);
var
I: Integer;
PtrD0,
PtrFX: PExtended;
begin
with
Temp, Option do
begin
//
Вычисление значений функции
PtrD0 := D0;
PtrFX :=
FX;
for
I := 0 to N do
begin
PtrFX^
:= Func(PtrD0);
PChar(PtrD0) := PChar(PtrD0) + N * SizeOf(Extended);
PChar(PtrFX) := PChar(PtrFX) + SizeOf(Extended);
end;
end;
end;
// Освежаем и сортируем FX + освежаем индексы
procedure
RefreshFX(var Temp: TNelderTempData; Option: TNelderOption);
var
I: Integer;
PtrIndexes:
PInteger;
begin
with
Temp, Option do
begin
//
Заполение индексов
PtrIndexes := Indexes;
for
I := 0 to N do
begin
PtrIndexes^ := I;
PChar(PtrIndexes) := PChar(PtrIndexes) + SizeOf(Integer);
end;
//
Сортировка
QuickSortFX(0, N, FX, Indexes);
//
Возвращаемое значение: Result := D0 + SizeOf(Extended) * Indexes[N]
PChar(PtrIndexes) := PChar(PtrIndexes) -
SizeOf(Integer);
XN :=
PExtended(PChar(D0) + N * SizeOf(Extended) * PtrIndexes^);
end;
end;
procedure
CalcC(var Temp: TNelderTempData; const Option: TNelderOption);
var
PtrC,
PtrD0: PExtended;
I, J: Integer;
begin
with
Temp, Option do
begin
PtrD0 :=
D0;
// C
:= 0;
FillChar(C^, N * SizeOf(Extended), 0);
// C
:= Sum (Xn)
for
I := 0 to N do
if
PtrD0 <> XN then
begin
PtrC
:= C;
for
J := 0 to N - 1 do
begin
PtrC^ := PtrC^ + PtrD0^;
PChar(PtrC) := PChar(PtrC) + SizeOf(Extended);
PChar(PtrD0) := PChar(PtrD0) + SizeOf(Extended);
end;
end
else
begin
//
Пропускаем вектор из D0:
PChar(PtrD0)
:= PChar(PtrD0) + N * SizeOf(Extended);
end;
// C
:= C / N
PtrC :=
C;
for
J := 0 to N - 1 do
begin
PtrC^
:= PtrC^ / N;
PChar(PtrC) := PChar(PtrC) + SizeOf(Extended);
end;
end;
end;
procedure
ReflectPoint(var Temp: TNelderTempData; const Option:
TNelderOption;
P:
PExtended; Factor: Extended);
var
PtrC,
PtrXN: PExtended;
I: Integer;
begin
with
Temp, Option do
begin
PtrXN :=
XN;
PtrC :=
C;
for
I := 0 to N - 1 do
begin
P^ :=
PtrC^ + Factor * (PtrC^ - PtrXN^);
PChar(P) := PChar(P) + SizeOf(Extended);
PChar(PtrC) := PChar(PtrC) + SizeOf(Extended);
PChar(PtrXN) := PChar(PtrXN) + SizeOf(Extended);
end;
end;
end;
const
SITUATION_EXPANSION = 0;
SITUATION_REFLECTION
= 1;
SITUATION_COMPRESSION = 2;
SITUATION_REDUCTION = 3;
function
DetectSituation(var Temp: TNelderTempData; const Option:
TNelderOption): Integer;
//FX, U:
PExtended; Func: PMathFunction; N: Integer; FU: Extended): Integer;
var
PtrFX:
PEXtended;
begin
with
Temp, Option do
begin
FU :=
Func(Temp.U);
if
FU < FX^ then
Result
:= SITUATION_EXPANSION
else
begin
PtrFX
:= PExtended(PChar(FX) + (N - 1) * SizeOf(Extended));
if
FU < PtrFX^ then
Result := SITUATION_REFLECTION
else
begin
PChar(PtrFX) := PChar(PtrFX) + SizeOf(Extended);
if
FU < PtrFX^ then
Result := SITUATION_COMPRESSION
else
Result := SITUATION_REDUCTION;
end;
end;
end;
end;
procedure
Expansion(var Temp: TNelderTempData; const Option:
TNelderOption);
var
FV:
EXtended;
LastIndex:
Integer;
TempPtr:
PChar;
begin
with
Temp, Option do
begin
ReflectPoint(Temp, Option, V, Gamma);
FV :=
Func(Temp.V);
// Коррекция FX
Move(FX^,
(PChar(FX) + SizeOf(Extended))^, N * SizeOf(Extended));
//
Заносим на первое место
if
FV < FU then
begin
FX^ :=
FV;
Move(V^, XN^, N * SizeOf(EXtended));
end
else
begin
FX^ :=
FU;
Move(U^, XN^, N * SizeOf(Extended));
end;
//
Коррекция Indexes
TempPtr
:= PChar(Indexes) + N * SizeOf(Integer);
LastIndex
:= PInteger(TempPtr)^;
Move(Indexes^, (PChar(Indexes) + SizeOf(Integer))^, N *
SizeOf(Integer));
Indexes^
:= LastIndex;
//
Коррекция XN
PChar(XN)
:= PChar(D0) + PInteger(TempPtr)^ * N * SizeOf(EXtended);
end;
end;
// Рекурсивный бинарный поиск точки, где будет
произведена вставка
// Интересно переделать в итерацию !!! (так
оптимальней)
function RecurseFind(FX:
PExtended; Value: Extended; L, R: Integer): Integer;
var
M: Integer;
Temp:
Extended;
begin
if R
< L then
begin
Result :=
L; // Result := -1 если поиск точный
Exit;
end;
M := (L +
R) div 2;
Temp :=
PExtended(PChar(FX) + M * SizeOf(Extended))^;
if
Temp = Value then
begin
Result :=
M;
Exit;
end;
if
Temp > Value then
Result :=
RecurseFind(FX, Value, L, M - 1)
else
Result :=
RecurseFind(FX, Value, M + 1, R)
end;
procedure
Reflection(var Temp: TNelderTempData; const Option:
TNelderOption);
var
InsertN:
Integer;
ShiftPosition: PChar;
//IndexesPtr:
PInteger;
//FV:
Extended;
//I:
Integer;
TempIndex:
Integer;
TempPtr:
PChar;
begin
with
Temp, Option do
begin
//
Определения позиции вставки
InsertN
:= RecurseFind(FX, FU, 0, N);
//
Сдвижка FX
ShiftPosition := PChar(FX) + InsertN * SizeOf(Extended);
Move(ShiftPosition^, (ShiftPosition + SizeOf(Extended))^,
(N -
InsertN) * SizeOf(Extended));
PExtended(ShiftPosition)^ := FU;
//
Коррекция D0
Move(U^,
XN^, N * SizeOf(EXtended));
//
Коррекция Indexes
TempPtr
:= PChar(Indexes) + N * SizeOf(Integer);
TempIndex
:= PInteger(TempPtr)^;
ShiftPosition := PChar(Indexes) + InsertN * SizeOf(Integer);
Move(ShiftPosition^, (ShiftPosition + SizeOf(Integer))^,
(N -
InsertN) * SizeOf(Integer));
PInteger(ShiftPosition)^ := TempIndex;
//
Коррекция XN
PChar(XN)
:= PChar(D0) + N * PInteger(TempPtr)^ * SizeOf(EXtended);
end;
end;
procedure
Compression(var Temp: TNelderTempData; const Option:
TNelderOption);
begin
with
Temp, Option do
begin
//
Вычисление новой точки
ReflectPoint(Temp, Option, U, Beta);
FU :=
Func(U);
Reflection(Temp, Option);
end;
end;
procedure
Reduction(var Temp: TNelderTempData; const Option:
TNelderOption);
var
ZeroPoint:
PExtended;
PtrD0,
PtrFX: PExtended;
PtrX0,
PtrX: PExtended;
FX0:
EXtended;
I, J:
Integer;
begin
with
Temp, Option do
begin
PChar(ZeroPoint) := PChar(D0) + N * Indexes^ * SizeOf(Extended);
PtrD0 :=
D0;
PtrFX :=
FX;
FX0 :=
FX^;
for
I := 0 to N do
begin
if
PtrD0 = ZeroPoint then
begin
//
Точка пропускается
PtrFX^ := FX0;
end
else
begin
//
Вычисляем точку:
PtrX
:= PtrD0;
PtrX0
:= ZeroPoint;
for
J := 0 to N - 1 do
begin
PtrX^ := 0.5 * (PtrX^ + PtrX0^);
PChar(PtrX)
:= PChar(PtrX) + SizeOf(Extended);
PChar(PtrX0) := PChar(PtrX0) + SizeOf(Extended);
end;
//
Вычисляем функцию
PtrFX^ := Func(PtrD0);
end;
PChar(PtrFX) := PChar(PtrFX) + SizeOf(Extended);
PChar(PtrD0)
:= PChar(PtrD0) + N * SizeOf(Extended);
end;
end;
RefreshFX(Temp, Option);
end;
var
It: Integer
= 0;
function
NeedStop(var Temp: TNelderTempData; const Option: TNelderOption):
Boolean;
var
PtrD0,
PtrC, PtrFX: PExtended;
Sum1, Sum2:
Extended;
I, J:
Integer;
begin
// Не
все параметры используются...
with
Temp, Option do
begin
Sum1 :=
0.0;
if
Delta > 0.0 then
begin
FC :=
Func(C);
PtrFX
:= FX;
for
I := 0 to N do
begin
Sum1
:= Sum1 + Sqr(PtrFX^ - FC);
PChar(PtrFX) := PChar(PtrFX) + SizeOf(Extended)
end;
Sum1 :=
Delta * Sqrt(Sum1 / (N + 1));
end;
Sum2 :=
0.0;
if
Delta < 1.0 then
begin
PtrD0
:= D0;
for
I := 0 to N do
begin
PtrC
:= C;
for
J := 0 to N - 1 do
begin
Sum2 := Sum2 + Sqr(PtrD0^ - PtrC^);
PChar(PtrC) := PChar(PtrC) + SizeOf(Extended);
PChar(PtrD0) := PChar(PtrD0) + SizeOf(Extended);
end;
end;
Sum2 := (1.0 - Delta) * Sqrt(Sum2 / (N +
1));
end;
Result :=
Sum1 + Sum2 < Eps;
end;
end;
procedure
Correct(var Temp: TNelderTempData; Option: TNelderOption);
var
S:
Extended;
PtrC:
PEXtended;
I: Integer;
begin
with
Temp, Option do
begin
S := (D1
+ (N - 1) * D2) / (N + 1);
PtrC :=
C;
for
I := 0 to N - 1 do
begin
PtrC^
:= PtrC^ - S;
PChar(PtrC) := PChar(PtrC) + SizeOf(Extended);
end;
BuildSimplex(Temp, Option);
MoveSimplex(Temp, Option, C);
end;
end;
function
Norm(X1, X2: PEXtended; N: Integer): Extended;
var
I: Integer;
begin
Result :=
0.0;
for
I := 0 to N - 1 do
begin
Result :=
Result + Sqr(X1^ - X2^);
PChar(X1)
:= PChar(X1) + SizeOf(Extended);
PChar(X2)
:= PChar(X2) + SizeOf(Extended);
end;
Result :=
Sqrt(Result);
end;
function
SolutionNelder(const Option: TNelderOption): Integer;
var
Temp:
TNelderTempData;
IsFirst:
Boolean;
begin
It := 0;
IsFirst :=
True;
Result :=
InitializeTempData(Temp, Option.N);
if
Result <> FIND_MIN_OK then
Exit;
try
// Шаг
№1: построение симплекса
BuildSimplex(Temp, Option);
// Шаг
№2: перенос симплекса в точку X0
MoveSimplex(Temp, Option, Option.X0);
repeat
//
Шаг №3: вычисление значений функции (+ сортировка)
CalcFX(Temp, Option);
RefreshFX(Temp, Option);
repeat
//
Шаг №4: вычисление центра тяжести без точки Indexes[N]
CalcC(Temp, Option);
//
Шаг №5: Вычисление точки U
ReflectPoint(Temp, Option, Temp.U,
Option.Alpha);
//
Шаг №6: Определение ситуации
Temp.FU := Option.Func(Temp.U);
case
DetectSituation(Temp, Option)
{Temp.FX,
Temp.U, Option.Func, Option.N, Temp.FU)} of
SITUATION_EXPANSION: // Растяжение
Expansion(Temp, Option);
SITUATION_REFLECTION:
Reflection(Temp, Option); // Отражение
SITUATION_COMPRESSION: // Сжатие
Compression(Temp, Option);
SITUATION_REDUCTION: // Редукция
Reduction(Temp, Option);
else
Assert(False, 'Другое не предусматривается');
end;
//
Шаг №7 критерий остановки
if
NeedStop(Temp, Option) then
Break;
until
False;
if not IsFirst then
begin
if
Norm(Option.X, Temp.C, Option.N) < Option.Eps then
begin
Move(Temp.C^, Option.X^, Option.N * SizeOf(Extended));
Break;
end;
end;
IsFirst
:= False;
Move(Temp.C^, Option.X^, Option.N * SizeOf(Extended));
case
Option.Mode of
FIND_MIN_MODE_STD: Break;
FIND_MIN_MODE_1: Correct(Temp, Option);
FIND_MIN_MODE_2:
begin
BuildSimplex(Temp, Option);
MoveSimplex(Temp, Option, Temp.C);
end;
end;
until
False;
Result :=
FIND_MIN_OK;
finally
FinalizeTempData(Temp);
end;
end;
function
FindMin_Nelder(const Option: TNelderOption): Integer;
var
UseOption:
TNelderOption;
begin
try
Result :=
CheckNelderOptionPtr(@Option);
if
Result <> FIND_MIN_OK then
Exit;
Result :=
CopyData(Option, UseOption);
if
Result <> FIND_MIN_OK then
Exit;
Result :=
SolutionNelder(UseOption);
finally
end;
end;
No comments:
Post a Comment