алгоритм Нелдера-Мида

[ciaonataha]

Мне надо использовать алгоритм Нелдера-Мида. Есть два варианта, какой из них правельнее?

http://www.kti.ru/data/84/m_n_mida.html

или

http://ru.wikipedia.org/wiki/%D0%9C%...B8%D0%B4%D0%B0

[ciaonataha]

Я взяла алгоритм из википедии и он не работает, не могу понять почему

//====================================================================================================
void simplex::amoeba(int nn,std::vector < std::vector<double> > &simplex,std::vector <double> &concentration, std::vector <double> &Y)
{
 int ihi, ilo, inhi, j;
 double fxr = 0.0, fxe = 0.0, fxs = 0.0;
 std::vector <double> midpoint(nn), next(nn), next_e(nn), next_s(nn);
 double fmin = 0.0;
 std::vector <double> fx(nn+1);
 evaluate_simplex(nn,simplex,concentration,fx);

 while (true)
  {
    simplex_extremes(nn,fx,&ihi,&ilo,&inhi);
    printf ("fx[ihi] = %e \n",fx[ihi]);
    simplex_bearings(nn,simplex,midpoint,ihi);
    if(check_tol(fx[ihi],fx[ilo],tol)) {printf ("zdes\n");break;}
    fxr = update_simplex(nn,ihi,simplex,midpoint,next,1.0,concentration);

      if (fxr < fx[ilo])
       {
        printf ("fxr < fx[ilo]\n");
        fxe = update_stretch(nn,next,midpoint,next_e,2.0,concentration);
          
          if (fxe < fx[ilo])
           {
            printf ("fxe < fx[ilo]\n"); 
            for (j = 0 ; j < nn; j++)
            {
              simplex[ihi][j] = next_e[j];
            }
           }
             
            if (fxe > fx[ilo])
             {
               printf ("fxe > fx[ilo]\n");
                for (j = 0; j < nn; j++)
                  {
                    simplex[ihi][j] = next[j];
                  }

             }
       }
 
     if ((fxr < fx[ihi])&&(fxr > fx[inhi]))
       {
         printf ("(fxr < fx[ihi])&&(fxr > fx[inhi])\n");
         for (j = 0; j < nn; j++)
          {
            double a = next[j];
            double b = simplex[ihi][j];
            next[j]  = b;
            simplex[ihi][j] = a;
          }     
            double a = fxr;
            double b = fx[ihi];
            fxr     = b;
            fx[ihi] = a; 

            fxs = update_simplex(nn,ihi,simplex,midpoint,next_s,-0.5,concentration);
       }
     if (fxr > fx[ihi])
       {
         printf ("fxr > fx[ihi]\n");
            fxs = update_simplex(nn,ihi,simplex,midpoint,next_s,-0.5,concentration);
            
       }

    if (fxs < fx[ihi])
       {
         printf ("fxs < fx[ihi]\n");
          for (j = 0; j < nn; j++)
           {
              simplex[ihi][j] = next_s[j];
           }
       }
          
    if (fxs > fx[ihi])
       {
         printf ("fxs > fx[ihi]\n");
         contract_simplex(nn,ilo,simplex,fx,concentration);
         //for (j = 0; j < nn+1; j++)
          //{printf ("%e ",fx[j]);}printf("\n");
       } 
    }
 
  //std::vector <double> Y(nn);
  for (j = 0; j < nn; j++)
   { 
     Y[j] = simplex[ilo][j];
     printf ("%e ",Y[j]);
   }printf ("\n");
 
}


}

функция:

(x*x+y-11.)*(x*x+y-11.)+(x+y*y-7.)*(x+y*y-7.)

[ciaonataha]

//====================================================================================================
void simplex::evaluate_simplex(int nn,std::vector < std::vector<double> > &simplex,std::vector <double> &concentration, std::vector <double> &fx)
{
 int i;
 std::vector <double> Y(nn);
 for (i = 0; i < nn+1; i++)
  { 
    Y[0]=simplex[i][0];Y[1]=simplex[i][1];Y[2]=simplex[i][2];Y[3]=simplex[i][3];Y[4]=simplex[i][4];Y[5]=simplex[i][5];Y[6]=simplex[i][6];Y[7]=simplex[i][7];
    fx[i] = (func->*nrfuncv)(nn,Y,concentration);
    printf ("%e %i\n",fx[i],i);
  }
}

//===================================================================================================
void simplex::simplex_extremes(int nn,std::vector <double> &fx,int *ihi, int *ilo, int *inhi)
{
 int i;

  if (fx[0] > fx[1])
   {
     //printf ("fx[0]>fx[1]\n");
     *ihi = 0; *ilo = *inhi = 1;
   }
   else
    {
      //printf ("fx[0]<fx[1]\n");
      *ihi = 1; *ilo = *inhi = 0;
    }
 
   for (i = 2; i < nn+1; i++)
    {
      if (fx[i]<=fx[*ilo])
       { 
         //printf ("fx[i]<fx[0]=%i\n",i);
         *ilo = i;
       }
        else if (fx[i] > fx[*ihi])
         {
           //printf ("fx[i]>fx[1]=%i\n",i);
           *inhi = *ihi;
           *ihi  = i;
         }
          else if (fx[i] > fx[*inhi])
            {
              //printf ("fx[i]>fx[0]=%i\n",i);
              *inhi = i;
            } 
    }    
 printf ("*ihi=%i *ilo=%i *inhi=%i \n",*ihi,*ilo,*inhi);
}

//=====================================================================================================
void simplex::simplex_bearings(int nn,std::vector < std::vector<double> > &simplex,std::vector <double> &midpoint,int ihi)
{
  int i, j;
  int alfa = 1;
  std::vector <double>  dummy(nn); 
  for (j = 0; j < nn; j++)
    {
      midpoint[j] = 0.0;
    }

 for (j = 0; j < nn; j++)
  {
    for (i = 0; i < nn+1; i++)
     {
      if (i != ihi)
       {
          dummy[j]     = simplex[i][j];
          midpoint[j] +=dummy[j];
       }
     }
  }
 midpoint.resize(nn);
  for (j = 0; j < nn; j++)
   {
     midpoint[j] /=nn;
   }
}

//======================================================================================================
int simplex::check_tol(double fmax, double fmin, double ftol)
{
 double delta = fabs(fmax-fmin);
 double accuracy = (fabs(fmax)+fabs(fmin))*ftol;
 return (delta < (accuracy + ZEPS));
}

//======================================================================================================
double simplex::update_simplex(int nn,int ihi,std::vector < std::vector<double> > &simplex,std::vector <double> &midpoint,std::vector <double> &next,double scale, std::vector <double> &concentration)
{
 int i;

 double fx;

 for (i = 0; i < nn; i++)
  {
   next[i] = (1+scale)*midpoint[i] - scale*simplex[ihi][i]; 
   printf ("next=%lf \n",next[i]);
  }
   fx   = (func->*nrfuncv)(nn,next,concentration);
   printf ("fxr = %lf \n",fx);
  
 return fx;
}

//==========================================================================================================================
double simplex::update_stretch(int nn,std::vector <double> &next,std::vector <double> &midpoint,std::vector <double> &next_e,double scale, std::vector <double> &concentration)
{
  int i;
  double fxe = 0.0;
 
  for (i = 0; i < nn; i++)
   {
      next_e[i] = (1-scale)*midpoint[i]+scale*next[i];
   }

   fxe = (func->*nrfuncv)(nn,next_e,concentration);
   printf ("fxe = %e \n",fxe);

 return fxe;
}
//===========================================================================================================================

void simplex::contract_simplex(int nn, int ilo,std::vector < std::vector<double> > &simplex,std::vector <double> &fx,std::vector <double> &concentration)
{
  int i, j;
  std::vector <double> Y(nn);  

  for (int i = 0; i < nn+1; i++)
    {
      if (i!=ilo)
       {
         for (j = 0; j < nn; j++)
          {
            simplex[i][j] = simplex[ilo][j]+ (simplex[i][j]-simplex[ilo][j])*0.5;
            //simplex[i][j] = (simplex[ilo][j]+simplex[i][j])*0.5;
            Y[0]=simplex[i][0];Y[1]=simplex[i][1];Y[2]=simplex[i][2];Y[3]=simplex[i][3];Y[4]=simplex[i][4];Y[5]=simplex[i][5];Y[6]=simplex[i][6];Y[7]=simplex[i][7];
    fx[i] = (func->*nrfuncv)(nn,Y,concentration);
          }
       }
    }