stogo/testfun.h

   1 #ifndef TESTFUN_H
   2 #define TESTFUN_H
   3
   4 #include "linalg.h"
   5 #include "tools.h"
   6 #include "stogo_config.h"
   7
   8 const double pi=fabs(acos(-1.));
   9
  10 /* The Matrix a and vector c are needed in the Shekel function */
  11 static double a[10][4]={ { 4,4,4,4 } ,
  12                          { 1,1,1,1 } ,
  13                          { 8,8,8,8 } ,
  14                          { 6,6,6,6 } ,
  15                          { 3,7,3,7 } ,
  16                          { 2,9,2,9 } ,
  17                          { 5,5,3,3 } ,
  18                          { 8,1,8,1 } ,
  19                          { 6,2,6,2 } ,
  20                          {7,3.6,7,3.6} };
  21 static double c[10]= { .1 , .2 , .2 , .4 , .4 , .6 , .3, .7 , .5 , .5 };
  22
  23 void Domain_Shekel(RTBox box) {
  24   box.lb=0.0 ; box.ub=10.0 ;
  25 }
  26
  27 double Objective_Shekel(RCRVector x) {
  28   int n=x.GetLength() ;
  29   double R=0.0, S;
  30   for(int i=0;i<10;i++) {
  31     S=0;
  32     for(int j=0;j<n;j++) S+=pow(x(j)-a[i][j],2);
  33     R-=1/(S+c[i]);
  34   }
  35   return R;
  36 }
  37
  38 void Gradient_Shekel(RCRVector x, RVector &grad) {
  39   int n=x.GetLength() ;
  40   double R;
  41   for(int k=0;k<n;k++) {
  42     R=0.0;
  43     for(int i=0;i<10;i++) {
  44       R+=(2.0*x(k)-2.0*a[i][k])/(pow(pow(x(0)-a[i][0],2.0)+pow(x(1)-a[i][1],2.0)
  45                                 +pow(x(2)-a[i][2],2.0)+pow(x(3)-a[i][3],2.0)+c[i],2.0));
  46     }
  47     grad(k)=R;
  48   }
  49 }
  50
  51
  52 /******************** Unimodal functions ******************/
  53
  54 void Domain_Rosenbrock(RTBox box) {
  55   box.lb=-10.0 ; box.ub=10.0 ;
  56 }
  57
  58 double Objective_Rosenbrock(RCRVector x) {
  59    double a=x(1)-x(0)*x(0) ;
  60    double b=1-x(0) ;
  61    return 100*a*a + b*b ;
  62 }
  63
  64 void Gradient_Rosenbrock(RCRVector x, RVector &grad) {
  65   grad(0)=200*(x(1)-x(0)*x(0))*(-2*x(0))-2*(1-x(0)) ;
  66   grad(1)=200*(x(1)-x(0)*x(0)) ;
  67 }
  68
  69
  70 void Domain_McCormick(RTBox box) {
  71   box.lb(0)=-1.5 ; box.ub(0)=4.0 ;
  72   box.lb(1)=-3.0 ; box.ub(1)=4.0 ;
  73 }
  74
  75 double Objective_McCormick(RCRVector x) {
  76   return sin(x(0)+x(1)) + pow(x(0)-x(1),2.0) - 1.5*x(0) + 2.5*x(1) + 1.0 ;
  77 }
  78
  79 void Gradient_McCormick(RCRVector x, RVector &grad) {
  80   grad(0)=cos(x(0)+x(1)) + 2*(x(0)-x(1)) - 1.5 ;
  81   grad(1)=cos(x(0)+x(1)) - 2*(x(0)-x(1)) + 2.5 ;
  82 }
  83
  84
  85 void Domain_BoxBetts(RTBox box) {
  86   box.lb(0)=0.9 ; box.ub(0)=1.2 ;
  87   box.lb(1)=9.0 ; box.ub(1)=11.2 ;
  88   box.lb(2)=0.9 ; box.ub(2)=1.2 ;
  89 }
  90
  91 double Objective_BoxBetts(RCRVector x) {
  92   double x0=x(0),x1=x(1),x2=x(2) ;
  93   double sum=0.0 ;
  94   for (int i=1 ; i<=10 ; i++)
  95     sum+=pow(exp(-0.1*i*x0)-exp(-0.1*i*x1)-(exp(-0.1*i)-exp(-1.0*i))*x2,2.0);
  96   return sum ;
  97 }
  98
  99 void Gradient_BoxBetts(RCRVector x, RVector &grad) {
 100   double x0=x(0),x1=x(1),x2=x(2) ;
 101   double g0=0.0, g1=0.0, g2=0.0 ;
 102   for (int i=1 ; i<=10 ; i++) {
 103     g0 += -0.2*(exp(-0.1*i*x0)-exp(-0.1*i*x1)
 104           -(exp(-0.1*i)-exp(-1.0*i))*x2)*i*exp(-0.1*i*x0);
 105     g1 += 0.2*(exp(-0.1*i*x0)-exp(-0.1*i*x1)-(exp(-0.1*i)
 106           -exp(-1.0*i))*x2)*i*exp(-0.1*i*x1);
 107     g2 += 2.0*(exp(-0.1*i*x0)-exp(-0.1*i*x1)
 108           -(exp(-0.1*i)-exp(-1.0*i))*x2)*(-exp(-0.1*i)+exp(-1.0*i));
 109   }
 110   grad(0)=g0 ; grad(1)=g1 ; grad(2)=g2 ;
 111 }
 112
 113
 114 void Domain_Paviani(RTBox box) {
 115  box.lb=2.001 ; box.ub=9.999 ;
 116 }
 117
 118 double Objective_Paviani(RCRVector x) {
 119   double a,b,sum=0.0, mul=1.0 ;
 120   int n=x.GetLength() ;
 121   for (int i=0 ; i<n ; i++) {
 122     a=log(x(i)-2.0) ; b=log(10.0-x(i)) ;
 123     sum+= a*a + b*b ;
 124     mul*= x(i) ;
 125   }
 126   return sum - pow(mul,0.2) ;
 127 }
 128
 129 void Gradient_Paviani(RCRVector x, RVector &grad) {
 130   double sum, mul=1.0 ;
 131   int n=10 ;
 132   for (int i=0 ; i<n ; i++) {
 133     mul*= x(i) ;
 134   }
 135
 136   for (int j=0 ; j<n ; j++) {
 137     sum=2*log(x(j)-2.0)/(x(j)-2.0) - 2*log(10.0-x(j))/(10.0-x(j)) ;
 138     grad(j) = sum - 0.2*(mul/x(j))/pow(mul,0.8) ;
 139   }
 140 }
 141
 142
 143 void Domain_Beale3(RTBox box) { // NB Make sure that there is only one minima
 144  box.lb=-0.5 ; box.ub=4.0 ;
 145 }
 146
 147 double Objective_Beale3(RCRVector x) {
 148   double x1=x(0), x2=x(1) ;
 149   double b1=1.500 + (x2 - 1.)*x1 ;
 150   double b2=2.250 + (x2*x2 - 1.)*x1 ;
 151   double b3=2.625 + (x2*x2*x2 - 1.)*x1 ;
 152   return b1*b1 + b2*b2 + b3*b3 ;
 153 }
 154
 155 void Gradient_Beale3(RCRVector x, RVector &grad) {
 156   double x1=x(0), x2=x(1) ;
 157   grad(0) = 2*(1.500 + (x2 - 1.)*x1)*(x2 - 1.)
 158           + 2*(2.250 + (x2*x2 - 1.)*x1)*(x2*x2 - 1.)
 159           + 2*(2.625 + (x2*x2*x2 - 1.)*x1)*(x2*x2*x2 - 1.) ;
 160
 161   grad(1) = 2*(1.500 + (x2 - 1.)*x1)*x1
 162           + 4*(2.250 + (x2*x2  - 1.)*x1)*x2*x1
 163           + 6*(2.625 + (x2*x2*x2 - 1.)*x1)*x2*x2*x1 ;
 164 }
 165
 166 /************** Difficult unimodal problems ****************/
 167 void Domain_Trid(RTBox box) {
 168   int n=(box.lb).GetLength() ;
 169   double tmp=pow(n,2.0);
 170   box.lb=-tmp ; box.ub=tmp ;  // [-n^2,n^2]
 171 }
 172
 173 double Objective_Trid(RCRVector x) {
 174   int n=x.GetLength();
 175   double sum1=0.0, sum2=0.0;
 176   for (int i=1 ; i<=n ; i++)
 177     sum1+=pow(x(i-1)-1,2.0);
 178   for (int i=2 ; i<=n ; i++)
 179     sum2+=x(i-1)*x(i-2);
 180   return sum1 - sum2;
 181 }
 182
 183 void Gradient_Trid(RCRVector x, RVector &grad) {
 184   int n=x.GetLength();
 185   grad(0)=2*(x(0)-1)-x(1) ;
 186   for (int i=1 ; i<=n-2 ; i++)
 187     grad(i)=2*(x(i)-1) - (x(i-1) + x(i+1)) ;
 188   grad(n-1)=2*(x(n-1)-1) - x(n-2) ;
 189 }
 190
 191 /******************** Multimodal functions ****************/
 192
 193 void Domain_OneDim(RTBox box) {
 194  box.lb=-25.0 ; box.ub=25.0 ;
 195 }
 196
 197 double Objective_OneDim(RCRVector x) {
 198   return x(0)*sin(x(0)) ;
 199 }
 200
 201 void Gradient_OneDim(RCRVector x, RVector &grad) {
 202   grad(0) = x(0)*cos(x(0)) + sin(x(0)) ;
 203 }
 204
 205 void Domain_Poly1(RTBox box) {
 206 // box.lb=-5.0 ; box.ub=5.0 ;
 207  box.lb=-4.0;box.ub=4.01;
 208 }
 209
 210 double Objective_Poly1(RCRVector x) {
 211   double a=(x(0)*x(0) + x(1)*x(1) - 11) ;
 212   double b=(x(0)+x(1)*x(1) - 7) ;
 213   return a*a + b*b ;
 214 }
 215
 216 void Gradient_Poly1(RCRVector x, RVector &grad) {
 217   double a=(x(0)*x(0) + x(1)*x(1) - 11) ;
 218   double b=(x(0)+x(1)*x(1) - 7) ;
 219   grad(0) = 4*x(0)* a + 2*b ;
 220   grad(1) = 4*x(1)* a + 4*x(1)*b ;
 221 }
 222
 223 void Domain_Jacobsen(RTBox box) {
 224  box.lb=-100.0 ; box.ub=100.0 ;
 225 }
 226
 227 double Objective_Jacobsen(RCRVector x) {
 228   double a=(x(0)*x(0) + x(1) - 11) ;
 229   double b=(x(0)+x(1)*x(1) - 7) ;
 230   return a*a + b*b + 3;
 231 }
 232
 233 void Gradient_Jacobsen(RCRVector x, RVector &grad) {
 234   double a=(x(0)*x(0) + x(1) - 11) ;
 235   double b=(x(0)+x(1)*x(1) - 7) ;
 236   grad(0) = 4*x(0)* a + 2*b ;
 237   grad(1) = 2*x(1)* a + 4*x(1)*b ;
 238 }
 239
 240 void Domain_Camel3(RTBox box) {
 241  box.lb=-1.0 ; box.ub=2.0 ;
 242 }
 243
 244 double Objective_Camel3(RCRVector x) {
 245   double a=x(0)*x(0) ;
 246   return 12*a - 6.3*a*a + a*a*a + 6*x(1)*(x(1)-x(0)) ;
 247 }
 248
 249 void Gradient_Camel3(RCRVector x, RVector &grad) {
 250   double a=x(0)*x(0) ;
 251   grad(0) = 24*x(0) - 25.2*a*x(0) + 6*a*a*x(0) - 6*x(1) ;
 252   grad(1) = 12*x(1) - 6*x(0) ;
 253 }
 254
 255
 256 void Domain_Hansen(RTBox box) {
 257   box.lb=-10.0 ; box.ub=10.0 ;
 258 }
 259
 260 double Objective_Hansen(RCRVector x) {
 261   return (cos(1.0)+2.0*cos(x(0)+2.0)+3.0*cos(2.0*x(0)+3.0)+4.0*cos(3.0*x(0)
 262         +4.0)+5.0*cos(4.0*x(0)+5.0))*(cos(2.0*x(1)+1.0)+2.0*cos(3.0*x(1)+2.0)
 263         +3.0*cos(4.0*x(1)+3.0)+4.0*cos(5.0*x(1)+4.0)+5.0*cos(6.0*x(1)+5.0));
 264 }
 265
 266 void Gradient_Hansen(RCRVector x, RVector &grad) {
 267    grad(0) = (-2.0*sin(x(0)+2.0)-6.0*sin(2.0*x(0)+3.0)-12.0*sin(3.0*x(0)+4.0)
 268            -20.0*sin(4.0*x(0)+5.0))*(cos(2.0*x(1)+1.0)+2.0*cos(3.0*x(1)+2.0)
 269            +3.0*cos(4.0*x(1)+3.0)+4.0*cos(5.0*x(1)+4.0)+5.0*cos(6.0*x(1)+5.0));
 270
 271    grad(1) = (cos(1.0)+2.0*cos(x(0)+2.0)+3.0*cos(2.0*x(0)+3.0)+4.0*cos(3.0*x(0)
 272              +4.0)+5.0*cos(4.0*x(0)+5.0))*(-2.0*sin(2.0*x(1)+1.0)
 273              -6.0*sin(3.0*x(1)+2.0)-12.0*sin(4.0*x(1)+3.0)
 274              -20.0*sin(5.0*x(1)+4.0)-30.0*sin(6.0*x(1)+5.0));
 275 }
 276
 277
 278 void Domain_Shubert(RTBox box) {
 279   box.lb=-10.0 ; box.ub=10.0 ;
 280 }
 281
 282 double Objective_Shubert(RCRVector x) {
 283    return -sin(2.0*x(0)+1.0)-2.0*sin(3.0*x(0)+2.0)-3.0*sin(4.0*x(0)+3.0)
 284           -4.0*sin(5.0*x(0)+4.0)-5.0*sin(6.0*x(0)+5.0)-sin(2.0*x(1)+1.0)
 285           -2.0*sin(3.0*x(1)+2.0)-3.0*sin(4.0*x(1)+3.0)-4.0*sin(5.0*x(1)+4.0)
 286           -5.0*sin(6.0*x(1)+5.0);
 287 }
 288
 289 void Gradient_Shubert(RCRVector x, RVector &grad) {
 290    grad(0) = -2.0*cos(2.0*x(0)+1.0)-6.0*cos(3.0*x(0)+2.0)-12.0*cos(4.0*x(0)+3.0)
 291              -20.0*cos(5.0*x(0)+4.0)-30.0*cos(6.0*x(0)+5.0);
 292    grad(1) = -2.0*cos(2.0*x(1)+1.0)-6.0*cos(3.0*x(1)+2.0)-12.0*cos(4.0*x(1)+3.0)
 293              -20.0*cos(5.0*x(1)+4.0)-30.0*cos(6.0*x(1)+5.0);
 294 }
 295
 296
 297 void Domain_Griewank(RTBox box) {
 298   box.lb=-500 ; box.ub=700 ;
 299 }
 300
 301 double Objective_Griewank(RCRVector x) {
 302   double sum=0 ;
 303   double prod=1 ;
 304   for (int i=0 ; i<10 ; i++) {
 305     sum+=x(i)*x(i) ;
 306     prod*=cos(x(i)/sqrt(double(i+1))) ;
 307   }
 308   return sum/4000.0-prod + 1 ;
 309 }
 310
 311 void Gradient_Griewank(RCRVector x, RVector &grad) {
 312   double prod=1 ;
 313   for (int i=0 ; i<10 ; i++) {
 314     prod*=cos(x(i)/sqrt(double(i+1))) ;
 315   }
 316   for (int i=0 ; i<10 ; i++) {
 317     grad(i)=x(i)/2000.0 + 1/sqrt(double(i+1))
 318            *prod/cos(x(i)/sqrt(double(i+1)))*sin(x(i)/sqrt(double(i+1))) ;
 319   }
 320 }
 321
 322
 323 void Domain_Levy(RTBox box) {
 324   int n=(box.lb).GetLength() ;
 325   switch (n) {
 326   case 4 :
 327     box.lb=-10.0 ; box.ub=10.0 ;
 328     break ;
 329   default:
 330     box.lb=-5.0 ; box.ub=5.0 ;
 331   }
 332 }
 333
 334 double Objective_Levy(RCRVector x) {
 335   int n=x.GetLength();
 336   double sum=0.0;
 337
 338   for (int i=0 ; i<=n-2 ; i++)
 339     sum+=pow(x(i)-1,2.0)*(1+pow(sin(3*pi*x(i+1)),2.0));
 340   return pow(sin(3*pi*x(0)),2.0) + sum + (x(n-1)-1)*(1+pow(sin(2*pi*x(n-1)),2.0));
 341 }
 342
 343
 344 void Gradient_Levy(RCRVector x, RVector &grad) {
 345   int n=x.GetLength();
 346
 347   grad(0)=6*sin(3*pi*x(0))*cos(3*pi*x(0))*pi + 2*(x(0)-1)*(1+pow(sin(3*pi*x(1)),2.0));
 348
 349   for (int i=1 ; i<=n-2 ; i++)
 350     grad(i)=6*pow(x(i-1)-1,2.0)*sin(3*pi*x(i))*cos(3*pi*x(i))*pi
 351       + 2*(x(i)-1)*(1+pow(sin(3*pi*x(i+1)),2.0)) ;
 352
 353   grad(n-1)=6*pow(x(n-2)-1,2.0)*sin(3*pi*x(n-1))*cos(3*pi*x(n-1))*pi
 354     + 1 + pow(sin(2*pi*x(n-1)),2.0) + 4*(x(n-1)-1)*sin(2*pi*x(n-1))*cos(2*pi*x(n-1))*pi;
 355 }
 356
 357
 358 void Domain_Kowalik(RTBox box) {
 359   box.lb=0 ; box.ub=5.0 ;
 360 }
 361
 362 double Objective_Kowalik(RCRVector x) {
 363   // First element in a and b is not used
 364   double a[]={999,0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246};
 365   double b[]={999,1./0.25,1./0.5,1.,1./2,1./4,1./6,1./8,1./10,1./12,1./14,1./16};
 366   double x1=x(0),x2=x(1),x3=x(2),x4=x(3);
 367
 368   double s=0;
 369   for (int i=1 ; i<=11 ; i++)
 370     s+=pow(a[i]- (x1*(b[i]*b[i] + b[i]*x2))/(b[i]*b[i] + b[i]*x3 + x4),2.0) ;
 371   return s;
 372 }
 373
 374 void Gradient_Kowalik(RCRVector x, RVector &grad) {
 375   double a[]={999,0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246};
 376   double b[]={999,1./0.25,1./0.5,1.,1./2,1./4,1./6,1./8,1./10,1./12,1./14,1./16};
 377   double x1=x(0),x2=x(1),x3=x(2),x4=x(3);
 378   double g1=0,g2=0,g3=0,g4=0,tmp;
 379   for (int i=1 ; i<=11 ; i++) {
 380     tmp=a[i]- (x1*(b[i]*b[i] + b[i]*x2))/(b[i]*b[i] + b[i]*x3 + x4) ;
 381     g1=g1-2*tmp*(b[i]*b[i]+b[i]*x2)/(b[i]*b[i]+b[i]*x3+x4) ;
 382     g2=g2-2*tmp*(x1*b[i])/(b[i]*b[i]+b[i]*x3+x4) ;
 383     g3=g3+2*tmp*x1*(b[i]*b[i]+b[i]*x2)*b[i]/pow(b[i]*b[i]+b[i]*x3+x4,2.0) ;
 384     g4=g4+2*tmp*x1*(b[i]*b[i]+b[i]*x2)/pow(b[i]*b[i]+b[i]*x3+x4,2.0) ;
 385   }
 386   grad(0)=g1; grad(1)=g2; grad(2)=g3; grad(3)=g4;
 387 }
 388
 389
 390 void Domain_Camel6(RTBox box) {
 391  box.lb=-5.0 ; box.ub=10.0 ;
 392 }
 393
 394 double Objective_Camel6(RCRVector x) {
 395   double x1=x(0),x2=x(1) ;
 396   return 4.0*x1*x1-0.21E1*pow(x1,4.0)+pow(x1,6.0)/3+x1*x2-4.0*x2*x2
 397     + 4.0*pow(x2,4.0);
 398 }
 399
 400 void Gradient_Camel6(RCRVector x, RVector &grad) {
 401   double x1=x(0),x2=x(1) ;
 402   grad(0) = 8.0*x1-0.84E1*x1*x1*x1+2.0*pow(x1,5.0)+x2;
 403   grad(1) = x1-8.0*x2+16.0*x2*x2*x2;
 404 }
 405
 406 /** Multimodal function with a huge number of local optima **/
 407
 408 void Domain_Rastrigin(RTBox box) {
 409   box.lb=-4.12 ; box.ub=6.12;
 410 }
 411
 412 double Objective_Rastrigin(RCRVector x) {
 413   double sum=0.0;
 414   int n=x.GetLength();
 415
 416   for (int i=1 ; i<=n ; i++)
 417      sum+=pow(x(i-1),2.0)-10*cos(2*pi*x(i-1))+10 ;
 418   return sum ;
 419 }
 420
 421 void Gradient_Rastrigin(RCRVector x, RVector &grad) {
 422   int n=x.GetLength();
 423
 424   for (int i=1 ; i<=n ; i++)
 425     grad(i-1)= 2*x(i-1) + 20*sin(2*pi*x(i-1))*pi ;
 426 }
 427
 428 void Domain_Schwefel(RTBox box) {
 429   box.lb=-500.0;box.ub=500.0;
 430 }
 431
 432 double Objective_Schwefel(RCRVector x) {
 433   double sum=0.0;
 434   int n=x.GetLength();
 435
 436   for (int i=0 ; i<n ; i++)
 437     sum+=x(i)*sin(sqrt(fabs(x(i))));
 438   return -sum;
 439 }
 440
 441 void Gradient_Schwefel(RCRVector x, RVector &grad) {
 442   int n=x.GetLength();
 443
 444   for (int i=0; i<n ; i++) {
 445     if (x(i)>=0)
 446       grad(i)=-( sin(sqrt(x(i)))+sqrt(x(i))/2.0*cos(sqrt(x(i))) );
 447     else
 448       grad(i)=-( sin(sqrt(-x(i)))+sqrt(-x(i))/2.0*cos(sqrt(-x(i))) );
 449   }
 450 }
 451
 452 /******* Difficult multimodal problem, PERM(n) *********/
 453
 454 void Domain_Perm(RTBox box) {
 455   int n=(box.lb).GetLength();
 456   box.lb=double(-n) ; box.ub=double(n) ;
 457 }
 458
 459 double ObjPerm(RCRVector x, double beta) {
 460   int n=x.GetLength();
 461   double s1=0.0;
 462   for (int k=1; k<=n ; k++) {
 463     double s2=0.0;
 464     for (int i=1 ; i<=n; i++)
 465       //  s2+=(i^k+beta)*((x(i-1)/i)^k-1);
 466       s2+=(pow(1.0*i,1.0*k)+beta)*(pow(x(i-1)/i,1.0*k)-1) ;
 467     s1+=s2*s2;
 468   }
 469   return s1;
 470 }
 471
 472 void GradPerm(RCRVector x, RVector &grad, double beta) {
 473   int n=x.GetLength();
 474   for (int j=1 ; j<=n ; j++) {
 475     double s1=0.0;
 476     for (int k=1 ; k<=n ; k++) {
 477       double s2=0.0;
 478       for (int i=1 ; i<=n; i++)
 479         //      s2+=(i^k+beta)*((x(i-1)/i)^k-1);
 480         s2+=(pow(1.0*i,1.0*k)+beta)*(pow(x(i-1)/i,1.0*k)-1);
 481       //    s1+=2*s2*(j^k+beta)/(x(j-1)*j^k)*k*x(j-1)^k;
 482       s1+=2*s2*(pow(1.0*j,1.0*k)+beta)/pow(1.0*j,1.0*k)*k*pow(x(j-1),k-1.0);
 483     }
 484     grad(j-1)=s1;
 485   }
 486 }
 487
 488 double Objective_Perm_4_50(RCRVector x) {
 489   return ObjPerm(x,50);
 490 }
 491
 492 void Gradient_Perm_4_50(RCRVector x, RVector &grad) {
 493   GradPerm(x,grad,50);
 494 }
 495
 496 double Objective_Perm_4_05(RCRVector x) {
 497   return ObjPerm(x,0.5);
 498 }
 499
 500 void Gradient_Perm_4_05(RCRVector x, RVector &grad) {
 501   GradPerm(x,grad,0.5);
 502 }
 503
 504 /******************* Powersum **************/
 505 void Domain_Powersum(RTBox box) {
 506   box.lb=0.0 ; box.ub=4.0 ;
 507 }
 508
 509 double Objective_Powersum(RCRVector x) {
 510   int n=x.GetLength();
 511   double b[]={8,18,44,114};
 512
 513   double s1=0.0;
 514   for (int k=1 ; k<=n ; k++) {
 515     double s2=0.0;
 516     for (int i=1 ; i<=n ; i++)
 517       s2+=pow(x(i-1),1.0*k);
 518     s1+=pow(s2-b[k-1],2.0);
 519   }
 520   return s1;
 521 }
 522
 523 void Gradient_Powersum(RCRVector x, RVector &grad) {
 524   int n=x.GetLength();
 525   double b[]={8,18,44,114};
 526
 527   for (int j=0 ; j<n ; j++) {
 528     double s1=0.0;
 529     for (int k=1 ; k<=n ; k++) {
 530       double s2=0.0;
 531       for (int i=1 ; i<=n ; i++)
 532         s2+=pow(x(i-1),1.0*k);
 533       s1+=2*(s2-b[k-1])*k*pow(x(j),k-1.0);
 534     }
 535     grad(j)=s1;
 536   }
 537 }
 538
 539 #endif