Оптимизация функции методом деформируемого многогранника (Метод Нелдера-Мида)

Previous  Top  Next

    
 

 

Code:

{ **** UBPFD *********** by kladovka.net.ru ****

>> Оптимизация функции методом деформируемого многогранника (Метод Нелдера-Мида)

 

Передаваемая структура:

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;

 

Возвращаемое значение - 0 если хорошо, иначе ошибка

 

Зависимости: Windows

Автор:       Mystic, mystic2000@newmail.ru, ICQ:125905046, Харьков

Copyright:   Mystic (посвящается Оксане в память о ее дипломе)

Дата:        25 апреля 2002 г.

********************************************** }

 

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;

 

©Drkb::04249