Search This Blog

Monday, November 4, 2013

Применение Метода градиентного спуска на С++, С#, Fotran, Delphi, Mathlab


Применение Метода градиентного спуска на С++, С#, Fotran, Delphi, Mathlab
Рассмотрим функцию f, считая для определенности, что она зависит от трех переменных x,y,z. Вычислим ее частные производные дf/дх, дf/ду, дf/дz и образуем с их помощью вектор, который называют градиентом функции:
grad f(x, у, z) = дf (х, у,z) /дх*i+дf( x, у, z)/ду*j+дf(x, y,z)/дг*k.
Здесь i, j, k - единичные векторы, параллельные координатным осям. Частные производные характеризуют изменение функции f по каждой независимой переменной в отдельности. Образованный с их помощью вектор градиента дает общее представление о поведении функции в окрестности точки (х, у,z). Направление этого вектора является направлением наиболее быстрого возрастания функции в данной точке. Противоположное ему направление, которое часто называют антиградиентным, представляет собой направление наиболее быстрого убывания функции. Модуль градиента 
                          ______________________                                  __ 
|grad (х, у,z)| =Ц (дf/дх (х, у,z))2 +(дf/ду( x, у, z))2+(дf/дг(x, y,z))2.
определяет скорость возрастания и убывания функции в направлении градиента и антиградиента. Для всех остальных направлений скорость изменения функции в точке (х, у, z) меньше модуля градиента. При переходе от одной точки к другой как направление градиента, так и его модуль, вообще говоря, меняются. Понятие градиента естественным образом переносится на функции любого числа переменных.
Перейдем к описанию метода градиентного спуска. Основная его идея состоит в том, чтобы двигаться к минимуму в направлении наиболее быстрого убывания функции, которое определяется антиградиентом. Эта идея реализуется следующим образом.
Выберем каким-либо способом начальную точку, вычислим в ней градиент рассматриваемой функции и сделаем небольшой шаг в обратном, антиградиентном направлении. В результате мы придем в точку, в которой значение функции будет меньше первоначального. В новой точке повторим процедуру: снова вычислим градиент функции и сделаем шаг в обратном направлении. Продолжая этот процесс, мы будем двигаться в сторону убывания функции. Специальный выбор направления движения на каждом шаге позволяет надеяться на то, что в данном случае приближение к наименьшему значению функции будет более быстрым, чем в методе покоординатного спуска.
Метод градиентного спуска требует вычисления градиента целевой функции на каждом шаге. Если она задана аналитически, то это, как правило, не проблема: для частных производных, определяющих градиент, можно получить явные формулы. В противном случае частные производные в нужных точках приходится вычислять приближенно, заменяя их соответствующими разностными отношениями:
дf ~  f(x1, ...,xi+Dxi, ..., xn) - f(x1, ..., xi, ..., xn)
 Отметим, что при таких расчетах Dxi ,нельзя брать слишком малым, а значения функции нужно вычислять с достаточно высокой степенью точности, иначе при вычислении разности
 Df(x1, ...,xi+Dxi, ..., xn) - f(x1, ..., xi, ..., xn)
 На рис. 3 изображены линии уровня той же функции двух переменных u= f (х, у), что и на рис. 1, и приведена траектория поиска ее минимума с помощью метода градиентного спуска.
Сравнение рис. 1 и 3 показывает, насколько более эффективным является метод градиентного спуска.


Posted by at No comments:

Метод покоординатного спуска С++,С#,Fotran,Delphi

Метод покоординатного спуска С++,С#,Fotran,Delphi

Пусть нужно найти наименьшее значение целевой функции u=f(M)=f(x1, x2, . . . ,xn). Здесь через М обозначена точка n-мерного пространства с координатами x1, x2, . . . ,xn: M=(x1, x2, . . . ,xn).Выберем какую-нибудь начальную точку М0=(x10, x20, . . . ,xn0) и рассмотрим функцию f при фиксированных значениях всех переменных, кроме первой: f(x1, x20,x30, . . . ,xn0 ). Тогда она превратится в функцию одной переменной  x1 . Изменяя эту переменную, будем двигаться от начальной точки x1=x10  в сторону убывания функции, пока не дойдем до ее минимума при x1=x11, после которого она начинает возрастать. Точку с координатами ( x11, x20,x30, . . . ,xn0) обозначим через М1, при этом f(M0) >= f(M1).

Фиксируем теперь переменные: x1=x11, x3= x30, . . . ,xn=xn0 и рассмотрим функцию f как функцию одной переменной  x2: f(x11, x22, x30 . . . ,xn0). Изменяя  x2 , будем опять двигаться от начального значения x2=x20   в сторону убывания функции, пока не дойдем до минимума при x2=x21 .Точку с координатами {x11, x21, x30 . . . xn0} обозначим через М2, при этом   f(M1) >=f(M2). 

Проведем такую же минимизацию целевой функции по переменным   x3, x4,  . . . ,xn. Дойдя до переменной xn, снова вернемся к x1 и продолжим процесс. Эта процедура вполне оправдывает название метода. С ее помощью мы построим последовательность точек М0,М1,М2, . . . , которой соответствует монотонная последовательность значений функцииf(M0) >= f (M1)>= f(M2) >=  ...    Обрывая ее на некотором шаге k можно приближенно принять  значение функции f(Mk) за ее наименьшее значение в рассматриваемой области.

Рис. 2. Поиск наименьшего значе-
ния функции методом покоординат-
ного спуска.

Отметим, что данный метод сводит задачу поиска наименьшего значения функции нескольких переменных к многократному решению одномерных задач оптимизации. Если целевая функция   f(x1, x2, ... ,xn)  задана явной формулой и является дифференцируемой, то мы можем вычисллить ее частные производные и использовать их для

определения направления убывания функции по каждой переменной и поиска соответствующих одномерных минимумов. В противном случае, когда явной формулы для целевой функции нет, одномерные задачи следует решать с помощью одномерных методов.

На рис. 2 изображены линии уровня некоторой функции двух переменных             u= f (х, у). Вдоль этих линий функция сохраняет постоянные значения, равные 1, 3, 5, 7, 9. Показана траектория поиска ее наименьшего значения, которое достигается в точке О, с помощью метода покоординатного спуска. При этом нужно ясно понимать, что рисунок служит только для иллюстрации метода. Когда мы приступаем к решению реальной задачи оптимизации, такого рисунка, содержащего в себе готовый ответ, у нас, конечно, нет.

 Пусть требуется решить задачу (2):

f(x) -->min,  х ОRn.                                      (2)


В двумерном пространтсве R2. Решение задачи (2) методом покоординатного спуска, иначе называемого методом Гаусса - Зейделя, производят по следующей общей схеме. 
Выбирают произвольно начальную точку х(0)   из   области определения функции f(х). Приближения х(k) определяются соотношениями (3): 
x(k+1)=x(k)+t(k)S(k)      (k=0,1,2, ...), 
где вектор направления спуска s(k)- это единичный вектор, совпадающий с каким-либо координатным направлением (например, если  S(k) параллелен х1, то S(k)= {1,0,0,...,0}, если он параллелен  x2, то S(k)={0, 1, 0, . . . ,0} и т.д.) ; величина t(k)   является решением задачи одномерной минимизации: 
f(x(k)+ts(k)) -->  min,   t ОR1,    (k=0,1,2, ...), 
и может определяться, в частности, методом сканирования. 
Детальная реализация общей схемы в двумерном случае R2 дает траекторий приближения к точке х* методом покоординатного спуска, состоящую из звеньев ломаной, соединяющих точки х(k), х~(k),x(k+1)   (k=0, 1,  2,)       (рис.2).    При k=0, исходя  из начальной точки х(0)=(x1(0),x2(0)), находят точку х~(0)= (x1~(0),x2(0)), минимума функции одной переменной f(x1,x2(0)); при этом f(x~(0))<=f(x(0)).Затем находят точку минимума x(1)  функции f (x1~(0),x2) по второй  координате. Далее делают следующий шаг вычислений при k=1. Полагают, что исходной точкой расчета является  х(1).   Фиксируя вторую  координату   точки  х(1),   находят   точку   минимума х~(1)= (x1~(1),x2(1)), функции  f(x1,x2(1)) одной переменной x(1); при  этом f(x~(1))<=f(x(1))<=f(x(0)). Точку х(2)   получают, минимизируя целевую функцию f(x1~(1),x2), вновь по коорданате х2, фиксируя координату x1~(1) ,точки x(1) , и т.д. 
Условием прекращения вычислительной процедуры при достижении заданной точности e может служить неравенство 
||x(k+1) - x(k) ||<e           (4) 
  
Блок-схема поиска минимума функции двух переменных методом покоординатного спуска. 


Posted by at No comments:

Метод деформируемого многогранника C++,C#,Delphi,Fotran (Метод Нелдера-Мида)


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;


Posted by at No comments:
Subscribe to: Posts (Atom)