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

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

{ **** 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, <a href="mailto:mystic2000@newmail.ru">mystic2000@newmail.ru</a>, 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;

Отправить комментарий

Проверка
Антиспам проверка
Image CAPTCHA
...