Модифицированный метод Хука-Дживса

[Voyageur]

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

[Dukales]

Не знаю, можно ли так делать по правилам форума, но вот исходники моей реализации метода Хука-Дживса для многомерного случая на C++ (в "C-стиле"): минимизация значения целевой функции многопараметрическим многомерным прямым поиском методом Хука-Дживса с ограничениями (и без них) с различным пределом шага для каждой параметра-константы. Если не нравятся "прямоугольные" ограничения, то можно использовать штрафные функции.
Сделаны по описанию алгоритма, взятому из книги Brian D.Bunday - "Basic optimisation methods". Код "многоканальный", если выражаться языком DSP. То есть можно серилизовано обращаться к данным из других потоков (расставлены критические секции) и возможно выполнять множество задач параллельно (для разных экземпляров структур HJD).
К ним ещё прилагался компилятор, который преобразует строку-выражение в код матсопроцессора (функция формата ENTRYPOINT в исходниках), но выложить его здесь - совсем оффтоп.
commonapi.h:

Используется синтаксис C++

#pragma once

#ifndef commonapiH
#define commonapiH

#include <assert.h>  
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <errno.h>
#include <ctype.h>
#include <malloc.h>
#include <time.h>
#include <algorith.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <windows.h>
#include <process.h>

#define ARRAYSIZE(ar) (sizeof(ar)/sizeof(ar[0]))

#endif

 

hook_and_jeeves.h:

Используется синтаксис C++

#ifndef hook_and_jeevesH
#define hook_and_jeevesH
//---------------------------------------------------------------------------

typedef __declspec(naked) double _fastcall (* ENTRYPOINT)(void * DSV, double DSC[]); // DSV: EAX, DSC: EDX, #pragma argsused

typedef struct _FUNC
{
  ENTRYPOINT EntryPoint;
  int ENTIRETYSV,
      ENTIRETYSC;
} FUNC, * pFUNC;

unsigned func_getvarcount(FUNC & func);

typedef struct _HJD_BASEPROPERTIES
{
  bool BordersAvailable;
  bool VariousPrecisions;
  volatile clock_t TimeOut;
  double Precision;
  double Divisor;
} HJD_BASEPROPERTIES, * pHJD_BASEPROPERTIES;

extern HJD_BASEPROPERTIES hjd_bp_defaults;

#pragma pack(push, 1)
typedef struct _BOUNDARY
{
  double Min,
         Max;
} BOUNDARY, * pBOUNDARY;
#pragma pack(pop)

#define HJD_SEARCH_ENABLED      0x00000000      // 000b
#define HJD_SEARCH_DISABLED     0x00000004      // 100b
#define HJD_SEARCH_STAND        0x00000000      // 000b
#define HJD_SEARCH_INC          0x00000001      // 001b
#define HJD_SEARCH_DEC          0x00000002      // 010b

typedef struct _HJD
{
  HJD_BASEPROPERTIES BP;
  volatile clock_t ExpendedTime;
  clock_t StartTime;
  ENTRYPOINT EntryPoint;
  size_t vsize;
  unsigned nco,
           pts;
  double Factor;
  double * BPV; // именно в BPV записан стартовый и конечный вектор
  double * UPV;
  void * Points;
  unsigned * MASK;
  double * Steps;
  double * StepsInitial;
  double * Precisions;
  BOUNDARY * Boundaries;
  long double OFV, BFV;
  unsigned FCC;
  CRITICAL_SECTION cs;
} HJD, * pHJD;

void hjd_init(HJD & hjd);
void hjd_free(HJD & hjd);

void hjd_setfunc(HJD & hjd, FUNC & f);

template <typename VARVECTOR> unsigned hjd_setpointsvec(HJD & hjd, VARVECTOR points[], unsigned pts)
{
  hjd.vsize = sizeof(VARVECTOR);

  hjd.pts = pts;

  unsigned size = pts * sizeof(VARVECTOR);
  hjd.Points = realloc(hjd.Points, size);
  for (unsigned i = 0; i < pts; i++) {
    ((VARVECTOR *)hjd.Points)[i] = points[i];
  }
  //memmove(hjd.Points, points, size);

  return (sizeof(VARVECTOR) / sizeof(double));
}

#define HJD_SET_STARTPOINT(_HJD, i, start) {(_HJD).BPV[(i)] = (start);}
#define HJD_SET_STARTVECTOR(_HJD, start) (memmove((_HJD).BPV, (start), (_HJD).nco * sizeof(double)))

#define HJD_SET_STEP(_HJD, i, step) {(_HJD).StepsInitial[(i)] = (step);}
#define HJD_SET_STEPS(_HJD, steps) (memmove((_HJD).StepsInitial, (steps), (_HJD).nco * sizeof(double)))

#define HJD_SET_PRECISION(_HJD, i, precision) {(_HJD).Precisions[(i)] = (precision);}
#define HJD_SET_PRECISIONS(_HJD, precisions) (memmove((_HJD).Precisions, (precisions), (_HJD).nco * sizeof(double)))

#define HJD_SET_BOUNDARY(_HJD, i, min, max) {(_HJD).Boundaries[(i)].Min = (min); (_HJD).Boundaries[(i)].Max = (max);}
#define HJD_SET_BOUNDARIES(_HJD, boundaries) (memmove((_HJD).Boundaries, (boundaries), (_HJD).nco * sizeof(BOUNDARY)))

#define HJD_MASK_BLOCK_POSITION(_HJD, position) {(_HJD).MASK[(position)] = HJD_SEARCH_DISABLED;}
#define HJD_MASK_UNBLOCK_POSITION(_HJD, position) {(_HJD).MASK[(position)] = HJD_SEARCH_ENABLED;}
#define HJD_SET_MASK_UNBLOCKALL(_HJD) {for (unsigned i = 0; i < (_HJD).nco; i++) (_HJD).MASK[i] = HJD_SEARCH_ENABLED;}

bool HookeAndJeevesTechnique(HJD & hjd);

//------------------------------------------------------------------------------
#endif

 

main.cpp:

Используется синтаксис C++

#include "commonapi.h"
#pragma hdrstop   
#include "hook_and_jeeves.h"

#pragma pack(push, 1)
typedef struct _DATAPOINT
{
  double _y,
         _x,
         _yc;
} DATAPOINT, * pDATAPOINT;
#pragma pack(pop)   

double _fastcall parabola(void * DSV, double DSC[])
{
  DATAPOINT & datapoint = *(DATAPOINT *)DSV;
  double temp = datapoint._x;
  datapoint._yc = DSC[0] + (DSC[1] + DSC[2] * temp) * temp;

  temp = datapoint._yc - datapoint._y;
  return temp * temp;
}

DATAPOINT Points[] = {{ 4.0, -2.0},
                      { 1.0, -1.0},
                      { 0.0,  0.0},
                      { 1.0,  1.0},
                      { 4.0,  2.0}};

const double BP[] = {0.0, 0.0, 0.0};
const double Steps[] = {0.001, 0.001, 0.001};
const double Precisions[] = {0.00001, 0.00001, 0.00001};
const BOUNDARY Boundaries[] = {{-10.0, 10.0},
                               {-10.0, 10.0},
                               {-10.0, 10.0}};

#pragma argsused
int main(int argc, char * argv[], char * env[])
{
  FUNC f;

  f.EntryPoint = parabola;
  f.ENTIRETYSV = 2;
  f.ENTIRETYSC = 2;
               
  HJD hjd;
  hjd_init(hjd);
  hjd_setfunc(hjd, f); 

  HJD_SET_STARTVECTOR(hjd, BP);
  HJD_SET_STEPS(hjd, Steps);
  HJD_SET_PRECISIONS(hjd, Precisions);
  HJD_SET_BOUNDARIES(hjd, Boundaries);

  HJD_SET_MASK_UNBLOCKALL(hjd);

  if (hjd_setpointsvec(hjd, Points, ARRAYSIZE(Points)) != func_getvarcount(f)) {
    throw;
  }

  HookeAndJeevesTechnique(hjd);

  for (unsigned i = 0; i < hjd.nco; i++) {
    printf("%lf\n", hjd.BPV[i]);
  }
  while (!kbhit()) continue;

  hjd_free(hjd);

  return EXIT_SUCCESS;
}

 

hook_and_jeeves.cpp:

Используется синтаксис C++

#include "commonapi.h"
#pragma hdrstop
#include "hook_and_jeeves.h"

unsigned func_getvarcount(FUNC & func)
{
  return (func.ENTIRETYSV + 1);
}

#define TWO_SECONDS 2000

HJD_BASEPROPERTIES hjd_bp_defaults =
{
  true,  // BordersAvailable
  true,  // VariousPrecisions
  LONG_MAX,  // TimeOut              // LONG_MAX ms // INFINITE
  0.0,   // Precision
  100.0 / 3.0   // Divisor
};

#pragma warn -rvl
#pragma argsused
static __declspec(naked) double _fastcall _retz(void * DSV, double DSC[])
{
  asm {
    FLDZ

    RET
  }
}
#pragma warn .rvl

void hjd_init(HJD & hjd)
{
  hjd.EntryPoint = _retz;

  hjd.nco = 0;
  hjd.Steps = (double *)calloc(0, sizeof(double));
  hjd.BPV = (double *)calloc(0, sizeof(double));
  hjd.UPV = (double *)calloc(0, sizeof(double));
  {
    hjd.StepsInitial = (double *)calloc(0, sizeof(double));
    hjd.Precisions = (double *)calloc(0, sizeof(double));
    hjd.Boundaries = (BOUNDARY *)calloc(0, sizeof(BOUNDARY));
  }
  hjd.MASK = (unsigned *)calloc(0, sizeof(unsigned));

  hjd.vsize = 0;
  hjd.pts = 0;
  hjd.Points = malloc(0);

  hjd.OFV = hjd.BFV = 0.0;
  hjd.ExpendedTime = hjd.StartTime = 0;
  hjd.BP = ::hjd_bp_defaults;

  InitializeCriticalSectionAndSpinCount(&hjd.cs, 0x400); //InitializeCriticalSection(&cs);
}

void hjd_setfunc(HJD & hjd, FUNC & f)
{
  hjd.EntryPoint = f.EntryPoint;

  hjd.nco = f.ENTIRETYSC + 1;
  unsigned size = hjd.nco * sizeof(double);
  hjd.Steps = (double *)realloc(hjd.Steps, size);
  hjd.BPV = (double *)realloc(hjd.BPV, size);
  hjd.UPV = (double *)realloc(hjd.UPV, size);
  {
    hjd.StepsInitial = (double *)realloc(hjd.StepsInitial, size);
    hjd.Precisions = (double *)realloc(hjd.Precisions, size);
    hjd.Boundaries = (BOUNDARY *)realloc(hjd.Boundaries, sizeof(BOUNDARY) * hjd.nco);
  }
  hjd.MASK = (unsigned *)realloc(hjd.MASK, sizeof(unsigned) * hjd.nco);
}

void hjd_free(HJD & hjd)
{
  DeleteCriticalSection(&hjd.cs);

  free(hjd.Points); hjd.Points = NULL;

  free(hjd.Steps); hjd.Steps = NULL;
  free(hjd.BPV); hjd.BPV = NULL;
  free(hjd.UPV); hjd.UPV = NULL;

  {
    free(hjd.StepsInitial); hjd.StepsInitial = NULL;
    free(hjd.Precisions); hjd.Precisions = NULL;
    free(hjd.Boundaries); hjd.Boundaries = NULL;
  }
  free(hjd.MASK); hjd.MASK = NULL;
}

static long double TargetFunction(HJD & hjd)
{
  long double TargetFunction = 0.0;

  for (unsigned i = 0; i < hjd.pts; i++) {
    TargetFunction += hjd.EntryPoint((unsigned __int8 *)hjd.Points + i * hjd.vsize, hjd.UPV);
  }

  hjd.FCC++;
  return TargetFunction;
}

static void SurveyRestricted(HJD & hjd)
{       
  for (unsigned i = 0; i < hjd.nco; i++) {
    if ((hjd.MASK[i] & HJD_SEARCH_DISABLED) == 0) {
      long double target_val;

      double displaced;
      double _ = hjd.UPV[i];

      //++
      displaced = _ + hjd.Steps[i];
      if (displaced < hjd.Boundaries[i].Max) {
        hjd.UPV[i] = displaced;

        target_val = TargetFunction(hjd);
        if (target_val < hjd.OFV) {
          hjd.OFV = target_val;
          hjd.MASK[i] = HJD_SEARCH_INC;

          continue;
        }
      }

      //--
      displaced = _ - hjd.Steps[i];
      if (displaced > hjd.Boundaries[i].Min) {
        hjd.UPV[i] = displaced;

        target_val = TargetFunction(hjd);
        if (target_val < hjd.OFV) {
          hjd.OFV = target_val;
          hjd.MASK[i] = HJD_SEARCH_DEC;

          continue;
        }
      }

      hjd.UPV[i] = _;
      hjd.MASK[i] = HJD_SEARCH_STAND;
    }
  }
}

static void Survey(HJD & hjd)
{
  for (unsigned i = 0; i < hjd.nco; i++) {
    if ((hjd.MASK[i] & HJD_SEARCH_DISABLED) == 0) {
      long double target_val;

      double _ = hjd.UPV[i];

      //++
      hjd.UPV[i] += hjd.Steps[i];

      target_val = TargetFunction(hjd);
      if (target_val < hjd.OFV) {
        hjd.OFV = target_val;
        hjd.MASK[i] = HJD_SEARCH_INC;

        continue;
      }

      //--
      hjd.UPV[i] = _ - hjd.Steps[i];

      target_val = TargetFunction(hjd);
      if (target_val < hjd.OFV) {
        hjd.OFV = target_val;
        hjd.MASK[i] = HJD_SEARCH_DEC;

        continue;
      }

      hjd.UPV[i] = _;
      hjd.MASK[i] = HJD_SEARCH_STAND;
    }
  }
}

static void PatternStepRestricted(HJD & hjd)
{
  for (unsigned i = 0; i < hjd.nco; i++) {
    switch (hjd.MASK[i]) {
      case HJD_SEARCH_INC : {
        hjd.UPV[i] += hjd.Steps[i];
        if (hjd.Boundaries[i].Max < hjd.UPV[i]) {
          hjd.UPV[i] = hjd.Boundaries[i].Max;
        }
      } break;
      case HJD_SEARCH_DEC : {
        hjd.UPV[i] -= hjd.Steps[i];
        if (hjd.Boundaries[i].Min > hjd.UPV[i]) {
          hjd.UPV[i] = hjd.Boundaries[i].Min;
        }
      }
    }
  }
}

static void PatternStep(HJD & hjd)
{
  for (unsigned i = 0; i < hjd.nco; i++) {
    switch (hjd.MASK[i]) {
      case HJD_SEARCH_INC : {
        hjd.UPV[i] += hjd.Steps[i];
      } break;
      case HJD_SEARCH_DEC : {
        hjd.UPV[i] -= hjd.Steps[i];
      }
    }
  }
}

static bool DecreaseStepsDistinguish(HJD & hjd)
{
  unsigned c = 0;

  for (unsigned i = 0; i < hjd.nco; i++) {
    hjd.Steps[i] *= hjd.Factor;
    if (hjd.Steps[i] <= hjd.Precisions[i]) {
      hjd.Steps[i] = hjd.Precisions[i];

      c++;
    }
  }

  return (c == hjd.nco);
}

static bool DecreaseStepsIndifference(HJD & hjd)
{
  unsigned c = 0;

  for (unsigned i = 0; i < hjd.nco; i++) {
    hjd.Steps[i] *= hjd.Factor;
    if (hjd.Steps[i] <= hjd.BP.Precision) {
      hjd.Steps[i] = hjd.BP.Precision;

      c++;
    }
  }

  return (c == hjd.nco);
}

bool HookeAndJeevesTechnique(HJD & hjd)
{
  bool hj = false;

  void (* _BEST_FIT_NEARBY)(HJD & hjd) = hjd.BP.BordersAvailable ? SurveyRestricted : Survey;
  void (* _PURSUE_BEST_DIRECTION)(HJD & hjd) = hjd.BP.BordersAvailable ? PatternStepRestricted : PatternStep;
  bool (* _DECREASE_STEPS_LENGTHS)(HJD & hjd) = hjd.BP.VariousPrecisions ? DecreaseStepsDistinguish : DecreaseStepsIndifference;

  EnterCriticalSection(&hjd.cs);

  hjd.StartTime = clock();
  hjd.ExpendedTime = 0;

  for (unsigned i = 0; i < hjd.nco; i++) {
    hjd.Steps[i] = hjd.StepsInitial[i];
  }

  hjd.Factor = 1.0 / hjd.BP.Divisor;

  hjd.FCC = 0;
  for (unsigned i = 0; i < hjd.nco; i++) {
    hjd.UPV[i] = hjd.BPV[i];
  }
  hjd.BFV = TargetFunction(hjd);

  do {
    hjd.OFV = hjd.BFV;
    for (unsigned i = 0; i < hjd.nco; i++) {
      hjd.UPV[i] = hjd.BPV[i];
    }

    _BEST_FIT_NEARBY(hjd);

    long double _ = hjd.OFV;

    if (hjd.BFV <= _) {
      if (hj) break;

      hj = _DECREASE_STEPS_LENGTHS(hjd);

      continue;
    }

    do {
      LeaveCriticalSection(&hjd.cs);
      hjd.ExpendedTime = clock() - hjd.StartTime;
      EnterCriticalSection(&hjd.cs);
      if (hjd.ExpendedTime > hjd.BP.TimeOut) {
        goto end;
      }

      hjd.BFV = _;
      for (unsigned i = 0; i < hjd.nco; i++) {
        hjd.BPV[i] = hjd.UPV[i];
      }

      _PURSUE_BEST_DIRECTION(hjd);

      _ = TargetFunction(hjd);
    } while (_ < hjd.BFV);
  } while (LeaveCriticalSection(&hjd.cs),
           hjd.ExpendedTime = clock() - hjd.StartTime,
           EnterCriticalSection(&hjd.cs),
           hjd.ExpendedTime < hjd.BP.TimeOut);

  end :

  hjd.ExpendedTime = clock() - hjd.StartTime;

  LeaveCriticalSection(&hjd.cs);

  return hj;
}

 

Так же есть альтернативная реализация на ассемблере (inline-asm Borland), но она работает раза в два медленнее.

-- Пн июн 13, 2011 14:21:15 --

Метод - именно модифицированный метод Хука-Дживса

[Dukales]

Небольшие изменения (давно писал то, что выше):
hooke_and_jeeves.h:

Используется синтаксис C++

#ifndef hooke_and_jeevesH
#define hooke_and_jeevesH
/*
  R. Hooke and T.A. Jeeves. Direct search solution of numerical and statistical
  problems. Journal of the Association for Computing Machinery (ACM), Volume 8 Issue 2, p. 212-229, April 1961
*/

typedef long double (* OBJECTIVE_FUNCTION)(double p[], void * param);

#pragma pack(1)
typedef struct _BOUNDARY
{
  double min;
  double max;
} BOUNDARY, * pBOUNDARY;
#pragma pack()

#define HJD_SEARCH_ENABLED      0x00000000
#define HJD_SEARCH_DISABLED     0x00000004
#define HJD_SEARCH_HOLD         0x00000000
#define HJD_SEARCH_INC          0x00000001
#define HJD_SEARCH_DEC          0x00000002

typedef struct _HJD
{
  bool R;                 // restriction on
  bool VP;                // various precisions on
  OBJECTIVE_FUNCTION of;  // Целевая функция
  void * param;           // Дополнительный параметр для целевой функции
  clock_t timeout;        // Таймаут
  double factor;          // Множитель для уменьшения значений шагов < 1
  unsigned nc;            // Количество параметров-констант целевой функции
  double * BPV;           // Именно в BPV записан стартовый и конечный вектор
  double * UPV;
  unsigned * mask;        // Маски: в каком направлении искать, блокировка поиска
  double * steps;         // Значения шагов. Можно после таймаута изменить
  double * precs;         // Предельные значения шагов. if VP
  double prec;            // Общее значение предельного значения шага для всех переменных. if !VP
  BOUNDARY * boundaries;  // Ограничения
  long double UOFV, BOFV; // Значение целевой функции в точке UPV и BPV
  unsigned OFCC;          // Количество вызовов целевой функции
} HJD, * pHJD;

void hjd_init(HJD & hjd); // Задаёт начальные значения параметров корректные для случайного вызова
void hjd_free(HJD & hjd); // Освобождает память

void hjd_setfunc(HJD & hjd, OBJECTIVE_FUNCTION of, void * param, unsigned nc, bool R, bool VP); // Назначает целевую функцию и вариант метода
void hjd_setperformance(HJD & hjd, clock_t timeout, double div, double prec = 0.0);             // Задаёт параметры, от которых зависит сходимость метода

#define HJD_SET_STARTVAR(_HJD, i, START) {(_HJD).BPV[(i)] = (START);}
#define HJD_SET_STARTVECTOR(_HJD, START) {memmove((_HJD).BPV, (START), (_HJD).nc * sizeof(double));}

#define HJD_SET_STEP(_HJD, i, STEP) {(_HJD).steps[(i)] = (STEP);}
#define HJD_SET_ALLSTEPS(_HJD, STEP) {for (unsigned i = 0; i < (_HJD).nc; i++) (_HJD).steps[i] = (STEP);} // Устанавливает значения всех шагов в одно и то же значение
#define HJD_IMPROVE_STEPS(_HJD, MULTIPLIER) {for (unsigned i = 0; i < (_HJD).nc; i++) (_HJD).steps[i] *= (MULTIPLIER);}
#define HJD_SET_STEPS(_HJD, STEPS) {memmove((_HJD).steps, (STEPS), (_HJD).nc * sizeof(double));}

#define HJD_SET_PRECISION(_HJD, i, PRECISION) {(_HJD).precs[(i)] = (PRECISION);}
#define HJD_SET_PRECISIONS(_HJD, PRECISIONS) {if (!(_HJD).VP) throw; memmove((_HJD).precs, (PRECISIONS), (_HJD).nc * sizeof(double));}

#define HJD_SET_BOUNDARY(_HJD, i, MIN, MAX) {(_HJD).boundaries[(i)].min = (MIN); (_HJD).boundaries[(i)].max = (MAX);}
#define HJD_SET_BOUNDARIES(_HJD, BOUNDARIES) {memmove((_HJD).boundaries, (BOUNDARIES), (_HJD).nc * sizeof(BOUNDARY));}

#define HJD_MASK_BLOCK_POSITION(_HJD, POSITION) {(_HJD).mask[(POSITION)] = HJD_SEARCH_DISABLED;} // Блокирует способность метода изменять параметр-константу
#define HJD_MASK_UNBLOCK_POSITION(_HJD, POSITION) {(_HJD).mask[(POSITION)] = HJD_SEARCH_ENABLED;}
#define HJD_SET_MASK_UNBLOCKALL(_HJD) {for (unsigned i = 0; i < (_HJD).nc; i++) (_HJD).mask[i] = HJD_SEARCH_ENABLED;}

clock_t HookeAndJeevesTechnique(HJD & hjd);

#endif
 

hooke_and_jeeves.cpp:

Используется синтаксис C++

#include "commonapi.h"
#pragma hdrstop
#include "hooke_and_jeeves.h"

#define objective_function(point, param) (hjd.OFCC++, hjd.of((point), (param)))

static bool BestNearby_R(HJD & hjd)
{
  hjd.UOFV = hjd.BOFV;
  for (unsigned i = 0; i < hjd.nc; i++) {
    hjd.UPV[i] = hjd.BPV[i];
  }

  for (unsigned i = 0; i < hjd.nc; i++) {
    if ((hjd.mask[i] & HJD_SEARCH_DISABLED) == 0) {
      long double OFV;

      double CCV = hjd.UPV[i];

      //++
      double D = CCV + hjd.steps[i];
      if (D < hjd.boundaries[i].max) {
        hjd.UPV[i] = D;
        OFV = objective_function(hjd.UPV, hjd.param);
        if (OFV < hjd.UOFV) {
          hjd.UOFV = OFV;
          hjd.mask[i] = HJD_SEARCH_INC;

          continue;
        }
      }

      //--
      D = CCV - hjd.steps[i];
      if (D > hjd.boundaries[i].min) {
        hjd.UPV[i] = D;
        OFV = objective_function(hjd.UPV, hjd.param);
        if (OFV < hjd.UOFV) {
          hjd.UOFV = OFV;
          hjd.mask[i] = HJD_SEARCH_DEC;

          continue;
        }
      }

      hjd.UPV[i] = CCV;
      hjd.mask[i] = HJD_SEARCH_HOLD;
    }
  }

  return (hjd.BOFV > hjd.UOFV);
}

static bool BestNearby_U(HJD & hjd)
{
  hjd.UOFV = hjd.BOFV;
  for (unsigned i = 0; i < hjd.nc; i++) {
    hjd.UPV[i] = hjd.BPV[i];
  }

  for (unsigned i = 0; i < hjd.nc; i++) {
    if ((hjd.mask[i] & HJD_SEARCH_DISABLED) == 0) {
      long double OFV;

      double CCV = hjd.UPV[i];

      //++
      hjd.UPV[i] += hjd.steps[i];
      OFV = objective_function(hjd.UPV, hjd.param);
      if (OFV < hjd.UOFV) {
        hjd.UOFV = OFV;
        hjd.mask[i] = HJD_SEARCH_INC;
      } else {
        //--
        hjd.UPV[i] = CCV - hjd.steps[i];
        OFV = objective_function(hjd.UPV, hjd.param);
        if (OFV < hjd.UOFV) {
          hjd.UOFV = OFV;
          hjd.mask[i] = HJD_SEARCH_DEC;
        } else {
          hjd.UPV[i] = CCV;
          hjd.mask[i] = HJD_SEARCH_HOLD;
        }
      }
    }
  }

  return (hjd.BOFV > hjd.UOFV);
}

static void PatternStep_R(HJD & hjd)
{
  hjd.BOFV = hjd.UOFV;
  for (unsigned i = 0; i < hjd.nc; i++) {
    hjd.BPV[i] = hjd.UPV[i];
  }

  for (unsigned i = 0; i < hjd.nc; i++) {
    switch (hjd.mask[i]) {
      case HJD_SEARCH_INC : {
        hjd.UPV[i] += hjd.steps[i];
        if (hjd.boundaries[i].max < hjd.UPV[i]) {
          hjd.UPV[i] = hjd.boundaries[i].max;
        }
      } break;
      case HJD_SEARCH_DEC : {
        hjd.UPV[i] -= hjd.steps[i];
        if (hjd.boundaries[i].min > hjd.UPV[i]) {
          hjd.UPV[i] = hjd.boundaries[i].min;
        }
      }
    }
  }

  hjd.UOFV = objective_function(hjd.UPV, hjd.param);
}

static void PatternStep_U(HJD & hjd)
{
  hjd.BOFV = hjd.UOFV;
  for (unsigned i = 0; i < hjd.nc; i++) {
    hjd.BPV[i] = hjd.UPV[i];
  }

  for (unsigned i = 0; i < hjd.nc; i++) {
    switch (hjd.mask[i]) {
      case HJD_SEARCH_INC : {
        hjd.UPV[i] += hjd.steps[i];
      } break;
      case HJD_SEARCH_DEC : {
        hjd.UPV[i] -= hjd.steps[i];
      }
    }
  }

  hjd.UOFV = objective_function(hjd.UPV, hjd.param);
}

static bool DecreaseSteps_V(HJD & hjd)
{
  unsigned nc = 0;

  for (unsigned i = 0; i < hjd.nc; i++) {
    hjd.steps[i] *= hjd.factor;
    if (hjd.steps[i] <= hjd.precs[i]) {
      hjd.steps[i] = hjd.precs[i];

      nc++;
    }
  }

  return (nc == hjd.nc);
}

static bool DecreaseSteps_G(HJD & hjd)
{
  unsigned nc = 0;

  for (unsigned i = 0; i < hjd.nc; i++) {
    hjd.steps[i] *= hjd.factor;
    if (hjd.steps[i] <= hjd.prec) {
      hjd.steps[i] = hjd.prec;

      nc++;
    }
  }

  return (nc == hjd.nc);
}

clock_t HookeAndJeevesTechnique(HJD & hjd)
{
  void (* PatternStep)(HJD & hjd) = hjd.R ? PatternStep_R : PatternStep_U;
  bool (* BestNearby)(HJD & hjd) = hjd.R ? BestNearby_R : BestNearby_U;
  bool (* DecreaseSteps)(HJD & hjd) = hjd.VP ? DecreaseSteps_V : DecreaseSteps_G;

  hjd.OFCC = 0;
  clock_t start = clock();

  hjd.UOFV = DBL_MAX;
  hjd.BOFV = objective_function(hjd.BPV, hjd.param);

  bool last = false;
  while (clock() - start < hjd.timeout) {
    if (hjd.BOFV > hjd.UOFV) {
      PatternStep(hjd);
    } else {
      while (!BestNearby(hjd)) {
        if (last) {
          return (clock() - start);
        } else {
          last = DecreaseSteps(hjd);
        }
      }
    }
  }

  return hjd.timeout;
}

#pragma warn -rvl
#pragma argsused
static __declspec(naked) long double _retz(double p[], void * param)
{
  asm {
    FLDZ

    RET
  }
}
#pragma warn .rvl

void hjd_init(HJD & hjd)
{
  hjd.of = _retz;

  hjd.nc = 0;
  hjd.BPV = NULL;
  hjd.UPV = NULL;
  {
    hjd.steps = NULL;
    hjd.precs = NULL;
    hjd.boundaries = NULL;
  }
  hjd.mask = NULL;
}

void hjd_setfunc(HJD & hjd, OBJECTIVE_FUNCTION of, void * param, unsigned nc, bool R, bool VP)
{
  hjd.R = R;
  hjd.VP = VP;

  hjd.of = of;
  hjd.param = param;

  hjd.nc = nc;
  unsigned size = nc * sizeof(double);
  hjd.BPV = (double *)realloc(hjd.BPV, size);
  hjd.UPV = (double *)realloc(hjd.UPV, size);
  {
    hjd.steps = (double *)realloc(hjd.steps, size);
    if (VP) {
      hjd.precs = (double *)realloc(hjd.precs, size);
    } else if (hjd.precs != NULL) {
      free(hjd.precs); hjd.precs = NULL;
    }
    if (R) {
      hjd.boundaries = (BOUNDARY *)realloc(hjd.boundaries, sizeof(BOUNDARY) * nc);
    } else if (hjd.boundaries != NULL) {
      free(hjd.boundaries); hjd.boundaries = NULL;
    }
  }
  hjd.mask = (unsigned *)realloc(hjd.mask, sizeof(unsigned) * nc);
}

void hjd_setperformance(HJD & hjd, clock_t timeout, double div, double prec)
{
  hjd.timeout = timeout;
  hjd.prec = prec;
  hjd.factor = 1.0 / div;
}

void hjd_free(HJD & hjd)
{
  free(hjd.BPV); hjd.BPV = NULL;
  free(hjd.UPV); hjd.UPV = NULL;
  {
    free(hjd.steps); hjd.steps = NULL;
    free(hjd.precs); hjd.precs = NULL;
    free(hjd.boundaries); hjd.boundaries = NULL;
  }
  free(hjd.mask); hjd.mask = NULL;
}

 

main.cpp:

Используется синтаксис C++

#include "commonapi.h"
#pragma hdrstop
#include "hooke_and_jeeves.h"

#pragma pack(1)
typedef struct _DATAPOINT
{
  double _y,
         _x,
         _yc;
} DATAPOINT, * pDATAPOINT;
#pragma pack()

long double _fastcall parabola(double c[], void * dp)
{
  DATAPOINT & datapoint = *(DATAPOINT *)dp;
  double temp = datapoint._x;
  datapoint._yc = c[2] + (c[1] + c[0] * temp) * temp;

  temp = datapoint._yc - datapoint._y;
  return temp * temp;
}

typedef struct _SEARCHSET
{
  void * points;
  unsigned pts;
  unsigned vsize;
} SEARCHSET, pSEARCHSET;

template <typename VARPOINT> inline void setvarvec(SEARCHSET & ss, VARPOINT * vv, unsigned pts)
{
  ss.vsize = sizeof(VARPOINT) / sizeof(double);
  ss.pts = pts;
  ss.points = vv;
}

static long double PARABOLA(double p[], void * param)
{
  SEARCHSET & ss = *(SEARCHSET *)param;

  long double P = 0.0;

  for (unsigned i = 0; i < ss.pts; i++) {
    P += parabola(p, (double *)ss.points + i * ss.vsize);
  }

  return P;
}

DATAPOINT points[] = {{ 5.0, -1.0},
                      { 2.0,  0.0},
                      { 1.0,  1.0},
                      { 2.0,  2.0},
                      { 5.0,  3.0}};

const double parabola_BP[] = {0.0, 0.0, 0.0};

const double precisions[] = {1E-6, 1E-6, 1E-5};
const BOUNDARY boundaries[] = {{-10.0, 10.0},
                               {-10.0, 10.0},
                               {-10.0, 10.0}};

/*
   Michael Powell,
   An Iterative Method for Finding Stationary Values of a Function
   of Several Variables,
   Computer Journal,
   Volume 5, 1962, pages 147-151.

   Powell's quartic function
   f(x) = (x1 + 10 * x2) ^ 2 + 5 * (x3 - x4) ^ 2 + (x2 - 2 * x3) ^ 4 + 10 * (x1 - x4) ^ 4
   start point: f(3, -1, 0, 1) = 215
   minimum: f(0, 0, 0, 0) = 0
*/
const double powell_BP[] = {3.0, -1.0, 0.0, 1.0};
#pragma argused
static long double powell(double p[], void * param)
{ 
  double a, b, c, d, t1, t2;

  a = p[0];
  b = p[1];
  c = p[2];
  d = p[3] ;

  t1  = a + 10.0 * b;
  t1 *= t1;
  t2  = c - d;
  t1 += 5.0 * t2 * t2;
  t2  = b - 2.0 * c;
  t2 *= t2;
  t1 += t2 * t2;
  t2  = a - d;
  t2 *= t2;

  return t1 + 10.0 * t2 * t2;
} 


// Rosenbrocks parabolic valley (banana function)
/*
   Rosenbrocks parabolic (banana) function
   f(x) = (1 - x1) ^ 2 + 100 * (x2 - x1 ^ 2) ^ 2
   start point: f(-1.2, 1) = 24.20
   minimum: f(1, 1) = 0
*/
const double rosen_BP[] = {-1.2, 1.0};
#pragma argused
static long double rosen(double p[], void * param)
{ 
  double a, b; 

  a = p[0];
  b = p[1] - a * a;
  a = 1.0 - a;
  return 100.0 * b * b + a * a;
}

// Schittowski problem 210: very difficult banana function
const double prob210_BP[] = {-1.2, 1.0};
#pragma argused
static long double prob210(double p[], void * param)
{
  double a, b; 

  a = p[0];
  b = p[1] - a * a;
  a = 1.0 - a;
  return 1.0E6 * b * b + a * a;
} 

// Schittowski problem 281: summation of ten cubics
const double prob281_BP[] = { 1.7528323e-02,
                              2.2503207e+00,
                              1.0325044e+00,
                              1.1047138e+00,
                              6.7813668e-01,
                              8.6220218e-01,
                              9.6336051e-01,
                              9.3734542e-01,
                              9.9586705e-01, 
                              9.5568268e-01};
#pragma argused
static long double prob281(double p[], void * param)
{ 
  double summ = 0.0; 
  for(unsigned i = 1; i <= 10; i++) 
  { 
    double temp = (p[i] - 1.0);
    summ += (i * i * i) * temp * temp;
  } 

  return pow(summ, 1.0 / 3.0);
}

#pragma argsused
int main(int argc, char * argv[], char * env[])
{
  printf("start hj\n");

  HJD hjd;
  hjd_init(hjd);

  #if false
  SEARCHSET ss;
  setvarvec(ss, points, ARRAYSIZE(points));

  hjd_setfunc(hjd, PARABOLA, &ss, ARRAYSIZE(parabola_BP), true, true);
  HJD_SET_PRECISIONS(hjd, precisions);
  HJD_SET_BOUNDARIES(hjd, boundaries);
  HJD_SET_MASK_UNBLOCKALL(hjd);

  hjd_setperformance(hjd, LONG_MAX, 56.399);
  HJD_SET_STARTVECTOR(hjd, parabola_BP);
  HJD_SET_ALLSTEPS(hjd, 1E-1);
  clock_t elapsed = HookeAndJeevesTechnique(hjd);
  printf("finish hj\n");
  printf("elapsed time = %d ms\n", elapsed);
  printf("OFCC = %d\n", hjd.OFCC);
  printf("OFV = %Lg\n", hjd.BOFV);

  printf("constants:\n");
  for (unsigned i = 0; i < hjd.nc; i++) {
    printf("\t% 3d: % 7lf +/- % 7lf\n", i, hjd.BPV[i], hjd.precs[i]);  
  }

  #elif true
  SEARCHSET ss;
  setvarvec(ss, points, ARRAYSIZE(points));

  hjd_setfunc(hjd, PARABOLA, &ss, ARRAYSIZE(parabola_BP), false, false);
  HJD_SET_MASK_UNBLOCKALL(hjd);

  FILE * f = fopen("divisor", "wb");
  for (double r = 2.0; r <= 100.0; r += 0.1) {
    hjd_setperformance(hjd, LONG_MAX, r, 0.0);
    HJD_SET_STARTVECTOR(hjd, parabola_BP);
    HJD_SET_ALLSTEPS(hjd, 1E-1);
    HookeAndJeevesTechnique(hjd);
    printf("%d\n", hjd.OFCC);
    fprintf(f, "%lg %d\n", r, hjd.OFCC);
  }
  fclose(f);

  putenv("PATH=C:\\Octave\\3.2.4_gcc-4.4.0\\bin");
  FILE * gnuplot = _popen("gnuplot.exe", "wb");
  fprintf(gnuplot, "%s\n", "plot 'divisor' using 1:2 with lines;");  
  fflush(gnuplot);
  while (!kbhit()) continue;
  fprintf(gnuplot, "%s\n", "exit;");
  fflush(gnuplot);
  _pclose(gnuplot);

  #elif false
  SEARCHSET ss;
  setvarvec(ss, points, ARRAYSIZE(points));

  hjd_setfunc(hjd, PARABOLA, &ss, ARRAYSIZE(parabola_BP), false, false);
  HJD_SET_MASK_UNBLOCKALL(hjd);

  FILE * f = fopen("precision", "wb");
  for (signed i = -1; i >= DBL_MIN_10_EXP; i--) {
    hjd_setperformance(hjd, LONG_MAX, 10.836, exp(M_LN10 * i));
    HJD_SET_STARTVECTOR(hjd, parabola_BP);       
    HJD_SET_ALLSTEPS(hjd, 1E-1);
    clock_t elapsed = HookeAndJeevesTechnique(hjd);
    //printf("elapsed time = %d ms\n", elapsed);
    printf("\tOFCC = %d\n", hjd.OFCC);
    //printf("OFV = %Lg\n", hjd.BOFV);
    //printf("general precision = %lg\n", hjd.prec);
    fprintf(f, "%i %d\n", -i, hjd.OFCC);
  } 
  fclose(f);  

  putenv("PATH=C:\\Octave\\3.2.4_gcc-4.4.0\\bin");
  FILE * gnuplot = _popen("gnuplot.exe", "wb");
  fprintf(gnuplot, "%s\n", "plot 'precision' using 1:2 with lines;");  
  fflush(gnuplot);
  while (!kbhit()) continue;
  fprintf(gnuplot, "%s\n", "exit;");
  fflush(gnuplot);
  _pclose(gnuplot);

  #elif false
  hjd_setfunc(hjd, powell, NULL, ARRAYSIZE(powell_BP), true, false);
  HJD_SET_BOUNDARIES(hjd, boundaries);
  HJD_SET_MASK_UNBLOCKALL(hjd);

  #define TIMEOUT 4000

  hjd_setperformance(hjd, TIMEOUT, 56.399, 1E-12);
  HJD_SET_STARTVECTOR(hjd, powell_BP);
  HJD_SET_ALLSTEPS(hjd, 1E-1);
  clock_t elapsed = HookeAndJeevesTechnique(hjd);
  while (elapsed >= TIMEOUT) {
    printf("OFCC = %d\n", hjd.OFCC);
    printf("OFV = %Lg\n", hjd.BOFV);  
    printf("general precision = %lg\n", hjd.prec);   
    printf("constants:\n");
    for (unsigned i = 0; i < hjd.nc; i++) {
      printf("\t% 3d: % 6lg\n", i, hjd.BPV[i]);
    }
    HJD_IMPROVE_STEPS(hjd, 2E2);
    printf("timeout, steps improved\n");
    elapsed = HookeAndJeevesTechnique(hjd);
  }
  printf("\tfinish!\n");
  if (TIMEOUT == LONG_MAX) {
    printf("elapsed time = %d ms\n", elapsed);
  }
  printf("OFCC = %d\n", hjd.OFCC);
  printf("OFV = %Lg\n", hjd.BOFV);  
  printf("general precision = %lg\n", hjd.prec);   
  printf("constants:\n");
  for (unsigned i = 0; i < hjd.nc; i++) {
    printf("\t% 3d: % 6lg\n", i, hjd.BPV[i]);
  }

  #elif false
  hjd_setfunc(hjd, rosen, NULL, ARRAYSIZE(rosen_BP), false, false); 
  HJD_SET_MASK_UNBLOCKALL(hjd);

  hjd_setperformance(hjd, LONG_MAX, 56.399, 0.0);
  HJD_SET_STARTVECTOR(hjd, rosen_BP);
  HJD_SET_ALLSTEPS(hjd, 1E-1);
  clock_t elapsed = HookeAndJeevesTechnique(hjd);
  printf("elapsed time = %d ms\n", elapsed);
  printf("OFCC = %d\n", hjd.OFCC);
  printf("OFV = %Lg\n", hjd.BOFV);

  printf("constants:\n");
  for (unsigned i = 0; i < hjd.nc; i++) {
    printf("\t% 3d: %lg +/- %lg\n", i, hjd.BPV[i], hjd.prec);  
  }   

  #elif false
  hjd_setfunc(hjd, prob210, NULL, ARRAYSIZE(prob210_BP), false, false);
  HJD_SET_MASK_UNBLOCKALL(hjd);

  hjd_setperformance(hjd, LONG_MAX, 56.399, 1.0E-9);
  HJD_SET_STARTVECTOR(hjd, prob210_BP);
  HJD_SET_ALLSTEPS(hjd, 1E-1);
  clock_t elapsed = HookeAndJeevesTechnique(hjd);
  printf("elapsed time = %d ms\n", elapsed);
  printf("OFCC = %d\n", hjd.OFCC);
  printf("OFV = %Lg\n", hjd.BOFV);

  printf("constants:\n");
  for (unsigned i = 0; i < hjd.nc; i++) {
    printf("\t% 3d: %lg +/- %lg\n", i, hjd.BPV[i], hjd.prec);  
  }

  #elif false
  hjd_setfunc(hjd, prob281, NULL, ARRAYSIZE(prob281_BP), false, false);
  HJD_SET_MASK_UNBLOCKALL(hjd);

  hjd_setperformance(hjd, LONG_MAX, 56.399, 1.0E-9);
  HJD_SET_STARTVECTOR(hjd, prob281_BP);
  HJD_SET_ALLSTEPS(hjd, 1E-1);
  clock_t elapsed = HookeAndJeevesTechnique(hjd);
  printf("elapsed time = %d ms\n", elapsed);
  printf("OFCC = %d\n", hjd.OFCC);
  printf("OFV = %Lg\n", hjd.BOFV);

  printf("constants:\n");
  for (unsigned i = 0; i < hjd.nc; i++) {
    printf("\t% 3d: %lg +/- %lg\n", i, hjd.BPV[i], hjd.prec);  
  }
  #endif

  while (!kbhit()) continue;

  hjd_free(hjd);

  return EXIT_SUCCESS;
}

 

commonapi.h:

Используется синтаксис C++

#pragma once

#ifndef commonapiH
#define commonapiH

#include <stdio.h>
#include <conio.h>
#include <stdlib.h>  
#include <mem.h>
#include <time.h>
#include <math.h>
#include <float.h>
#include <limits.h>

#define ARRAYSIZE(ar) (sizeof(ar)/sizeof(ar[0]))

#endif