snobfit.snobfit

50 """ 51 Minimization of a function over a box in R^n 52 53 (n = dimension of the problem) 54 """ 55 56 # The goldenSection number ( 1/2*[numpy.sqrt(5)-1] ) 57 goldenSection = 0.618 58

59 - def __init__(self, 60 func, 61 x0, 62 bounds, 63 dx=None, 64 constraint=None, 65 p=0.8, 66 dn=5, 67 fglob=None, 68 maxiter=1000, 69 maxfun=1000, 70 disp=0, 71 retall=0, 72 seed=None, 73 xglob=None, 74 nstop=250, 75 fac=0.0, 76 xtol=1.e-6, 77 ftol=1.e-6, 78 rtol=1.e-6, 79 isLeastSqrt = False, 80 retuct=False, 81 callback=None 82 ):

83 """ 84 Inputs: 85 ------ 86 For properly use Snobfit, the user must input the follow variables: 87 func the Python function or method to be minimized. 88 x0 the initial guess. 89 bounds the box boundary, it is a list of (lowBounds, highBounds). 90 91 Additional Inputs: 92 ------------------ 93 p probability of generating a evaluation point of Class 4. 94 fglob the user specified global function value. 95 xglob the user specified global minimum. 96 nstop number of times no improvement is tolerated 97 rtol a relative error for checking stopping criterion. 98 maxiter the maximum number of iterations to perform. 99 disp non-zero if fval and warnflag outputs are desired. 100 retall non-zero to return list of solutions at each iteration. 101 102 dx resolution vector, i.e. the i_th coordinate of a point to be 103 generated is an integer-valued multiple of dx[i] 104 105 Only used for the definition of a new problem n-vector of 106 minimal steps, i.e., two points are considered to be different 107 if they differ by at least dx[i] in coordinate i 108 109 dn In the presence of noise, fitting reliable linear models near 110 some point requires the use of a few, say dn more data points 111 than parameters. ( default 5 ) 112 113 seed The seed for numpy.random. If the user set this number, it makes 114 sure that it returns the same solution from repeated tests 115 116 xtol acceptable relative error in xopt for convergence. 117 ftol acceptable relative error in func(xopt) for convergence. 118 maxfun the maximum number of function evaluations. 119 120 isLeastSqrt the minimisation uses the least-squares function or not 121 retuct return uncertainty or not 122 123 callback an optional user-supplied function to call after each 124 iteration. It is called as callback(n,xbest,fbest,improved) 125 126 constraint an instance of SoftConstraints 127 128 Note: At this stage, we don't implement the code use xtol, ftol, 129 and maxfun. Later we should use it. 130 """ 131 # The initial guess of the fitting 132 self.x0 = x0 133 134 # The dimension of the problem 135 self.n = len(self.x0) 136 137 # User defined number 138 self.dn = dn 139 140 # The suggested number of the evaluation points for fitting. 141 self.nreq = self.n + self.dn + 1 142 143 # Number of safeguarded nearest neighbors 144 self.snn = self.n + self.dn 145 146 # the Python function or method to be minimized. 147 self.func = func 148 149 # the constraint 150 self.constraint = constraint 151 152 # The user defined lower and upper box bounds 153 (self.u1, self.v1) = bounds 154 155 # resolution vector 156 self.dx = dx 157 158 # Probability of generating recommended evaluation points of Class 4 159 self.p = p 160 161 # non-zero if fval and warnflag outputs are desired. 162 self.disp = disp 163 164 # non-zero to return list of solutions at each iteration 165 self.retall = retall 166 167 # set the global values 168 self.xglob = xglob 169 self.fglob = fglob 170 171 # check stopping criterion (best function value is found 172 # with a relative error < rtol 173 self.rtol = rtol 174 self.ftol = ftol 175 self.xtol = xtol 176 177 # number of times no improvement is tolerated 178 self.nstop = nstop 179 180 # the factor for multiplicative perturbation of the data 181 self.fac = fac 182 183 # the maximum number of iterations to perform. 184 self.maxiter = maxiter 185 186 # the maximum number of function evaluations. 187 self.maxfun = maxfun 188 189 # user-supplied function to call after each iteration. 190 self.callback = callback 191 192 # If the minimisation uses the least-squares function 193 # eg. f_i = (Y(x, t_i) - y_i) / \sigma_i or not 194 self.isLeastSqrt = isLeastSqrt 195 196 # Return the uncertainty the fitting parameters or not? 197 self.retuct = retuct 198 199 # Set the seed for numpy.random 200 self.seed = seed 201 if self.seed is not None: 202 numpy.random.seed( [self.seed] ) 203 204 # The rows are a set of new points entering the optimization 205 # algorithm together with their function values 206 self.x = self._setInitRecommendedEvalPoints( x0 ) 207 208 # The corresponding function values ( for self.x ) 209 # and their uncertainties. 210 # a value df[i] <= 0 indicates that the corresponding uncertainty 211 # is not known, and the program resets it to numpy.sqrt(eps) 212 (self.f, self.df) = self._setInitial_fdf() 213 214 # Current number of the recommended evaluation points 215 self.nx = self.x.shape[0] 216 217 # Other initialization 218 self.setInit()

219 220

221 - def _setInitRecommendedEvalPoints(self, x0):

222 """ 223 From the initial guess, we set up the first set of 224 the recommended evaluation points 225 226 Note: Make sure the initial guess x0 is 227 in recommended evaluation points. 228 """ 229 x = numpy.dot( rand(self.nreq-1, self.n), 230 numpy.diag(self.v1-self.u1) 231 ) + \ 232 numpy.array( [self.u1] ).repeat(self.nreq-1,axis=0) 233 234 x = numpy.append( [x0], x, axis=0 ) 235 return x

236 237

238 - def _setInitial_fdf(self):

239 """ 240 Computation of the function values of the first set of the recommended 241 evaluation function (if necessary, with additive noise) 242 """ 243 f = vector(self.nreq) 244 df = vector(self.nreq) 245 if self.constraint is not None: 246 self.F1 = self.constraint.F1() 247 self.F2 = self.constraint.F2() 248 fm = vector(self.nreq) 249 f0 = numpy.inf 250 251 for i in xrange(self.nreq): 252 f[i] = self.func( self.x[i,:] ) 253 df[i] = max(3*self.fac, numpy.sqrt(eps)) 254 255 if self.constraint is not None: 256 FF = self.constraint.F(self.x[i,:]) 257 if (sum( numpy.logical_and( self.F1<= FF, 258 FF <= self.F2) ) == self.n ): 259 f0 = min(f0, f[i]) 260 261 262 if self.constraint is not None: 263 Delta = numpy.median( abs(f-f0) ) 264 for i in xrange(self.nreq): 265 fm[i] = self.softmerit(f[i], f0, Delta, self.x[i,:]) 266 267 # Save it for later use 268 self._f0 = f0 269 self._Delta = Delta 270 return (fm,df) 271 272 return (f,df)

273 274

275 - def setInit(self):

276 """ 277 Some other initialization of Snobfit. 278 """ 279 if len(self.f)==0: 280 self.f = vector() 281 self.df = vector() 282 283 # If the uncertainty of the function value is negative, 284 # set it be the square root of the machine precision. 285 ind = find( self.df <= 0 ) 286 if len(ind) != 0: 287 self.df[ind] = numpy.sqrt(eps) 288 289 # Defines the vector of minimal distances between two 290 # points suggested in a single call to Snobfit 291 self.dy = 0.1*(self.v1-self.u1) 292 293 if len( find(self.dx==0) ) != 0 : 294 raise ValueError('dx should only contain positive entries') 295 296 # ---------------------------------------------------------- 297 # Some other common data block 298 # ---------------------------------------------------------- 299 # The search bounds [u,v] 300 self.u = self.u1.copy() 301 self.v = self.v1.copy() 302 303 # vector containing the pointers to the points 304 # where the function value could not be obtained 305 self.fnan = numpy.array([]) 306 307 # near[i,:] is a vector containing the indices of 308 # the nearest neighbors of the point x[i,:] 309 self.near = self.near = numpy.zeros( (self.nx,self.snn), 'int' ) 310 311 # d[j] is maximal distance between x[j,:] and one of its neighbors, 312 # or Initial value for the diameter of the trust region 313 # for local minimization from the best point 314 self.d = inf * numpy.ones(self.nx) 315 316 # small[j] is an logarithmic measure of box j 317 self.small = ivector() 318 319 # nsplit[j,i] number of times box j has been split along ith coordinate 320 self.nsplit = ivector() 321 322 # nxFval[j] is the number of times the function value 323 # of x[j,:] has been measured 324 self.nxFval = None 325 326 # t[j] is nxFval[j] times the variance of the function 327 # values measured for the point x[j,:] 328 self.t = None 329 330 # At the end of fit, ther should be nreq x (n+2)-matrix 331 self.request = numpy.zeros( (0, self.n+2) ) 332 333 # Initialize J4, which contains the indices of boxes marked for 334 # the generation of the recommended evaluation points point of Class 4, 335 # as the empty set. 336 self.J4 = ivector() 337 338 # Pointer pointing to the new boxes and 339 # boxes whose nearest neighbors have changed 340 self.inew = ivector() 341 342 # y, fy are a potential evaluation point and its function values 343 self.y = numpy.zeros( (self.nx, self.n) ) 344 self.fy = vector( self.nx ) 345 346 # The gradient of model 347 self.gradient = numpy.zeros( (self.nx, self.n) ) 348 349 # The standard deviation of model errors 350 self.sigma = vector( self.nx ) 351 352 # The follow variables are used to improve the speed. 353 self.dx2 = self.dx**2

354 355

356 - def _warning():

357 """ Warning message from snobfit """ 358 raise ValueError, """WARNING: The algorithm was not able to generate \ 359 the desired number of points.\n 360 Change the search region or refine resolution vector dx""" 361 sys.exit(1)

362 363

364 - def _calc_DMat(self,f):

365 """ Compute the D Matrix """ 366 return numpy.diag( f /self.dx2 ) # self.dx2 = self.dx**2

367 368

369 - def _calc_Smallness(self, xu, xl, v, u ):

370 """ Compute the smallness 371 372 This quantity roughly measures how many bisections are 373 necessary to obtain this box from [u, v]. 374 We have S = 0 for the exploration box [u, v], 375 and S is large for small boxes. 376 For a more thorough global search, boxes with low smallness 377 should be preferably explored when selecting new points for evaluation. 378 """ 379 small=-numpy.sum(numpy.round(numpy.log2((xu-xl)/ 380 numpy.array([v-u]).repeat(xu.shape[0],axis=0))), 381 axis=1 382 ) 383 return toInt( small )

384 385

386 - def _addRequest(self, newReq):

387 """ Add a new request. """ 388 self.request = numpy.append( self.request, newReq, axis=0 )

389 390

391 - def _isMarkedForClass4( self, i, xu, xl, v, u ):

392 """ 393 Decide the index of box ( i.e., i) is marked for the generation of 394 recommended evaluation points of Class 4 or not 395 """ 396 x = (xu[i,:]-xl[i,:])/(v-u) 397 if numpy.min(x) <= 0.05*numpy.min(x): 398 return True 399 else: 400 return False

401 402

403 - def _sumFind(self, x, n=0, axis=1):

404 """ A helper function """ 405 return find( numpy.sum(x, axis=axis) == n )

406 407

408 - def _maxmin(self, v):

409 """ A helper function """ 410 return ( max(v), min(v) )

411 412

413 - def _estimate(self, f0, df0, y, x, g, sigma):

414 """ 415 Calculate the Formula (2) of SNOBFIT by Arnold Neumaier 416 417 f0, df0: function value and its uncertainty in point y 418 """ 419 z = y-x 420 D = self._calc_DMat( df0 ) 421 f = f0 + numpy.dot(z, g) + sigma*( dot3(z, D, z.T) + df0 ) 422 return float(f)

423 424

425 - def _quadEstimate(self, f0, y, x, g, G):

426 """ 427 Calculate the Formula (3) of SNOBFIT by Arnold Neumaier 428 """ 429 z = y-x 430 f = f0 + numpy.dot(z, g) + 0.5*dot3(z, G, z.T) 431 return float(f)

432 433

434 - def newRequest(self, x, f, i):

435 _req = numpy.array([]) 436 _req = numpy.append( _req, x) 437 _req = numpy.append( _req, f) 438 _req = numpy.append( _req, i) 439 return _req

440 441

442 - def round(self, 443 x, # vector of length n 444 u, # lower bounds of x 445 v # upper bounds of x 446 ):

447 """ 448 A point x is projected into the interior of [u,v] and x[i] is 449 rounded to the nearest integer multiple of self.dx[i] 450 451 Return the projected and rounded version of x 452 """ 453 x = numpy.minimum(numpy.maximum(x,u),v) 454 x = numpy.round(x/self.dx)*self.dx 455 456 i1 = find(x<u) 457 i2 = find(x>v) 458 if len(i1) != 0: x[i1] += self.dx[i1] 459 if len(i2) != 0: x[i2] -= self.dx[i2] 460 461 return x

462 463

464 - def rsort(self, x):

465 """ 466 Sort x in increasing order, remove multiple entries. 467 468 Warning: when you use this function, make sure x and w is row vector 469 """ 470 x = numpy.sort(x) 471 472 # Remove repetitions 473 n = len(x) 474 xnew = numpy.append( x[1:n], inf ) 475 ind = find( xnew != x ) 476 x=x[ind] 477 478 return x

479 480 481 # --------------------------------------------------

482 - def quadMin(self, 483 a, b, # coefficients of the quadratic polynomial 484 xl, xu # bounds (xl < xu) 485 ):

486 """ 487 Minimization of the quadratic polynomial 488 p(x) = a*x^2+b*x over [xl,xu] 489 490 Output: 491 x minimizer in [xl,xu] 492 """ 493 if a > 0: 494 x = -0.5*b/a 495 x = numpy.minimum( numpy.maximum(xl,x), xu ) 496 497 else: 498 fl = a*xl**2 + b*xl 499 fu = a*xu**2 + b*xu 500 if fu <= fl: 501 x = xu 502 else: 503 x = xl 504 505 return x

506 507

508 - def _newBounds(self, 509 xnew, # the rows of xnew contain the new points 510 xl, # xl[j,:] is the lower bound of box j 511 xu, # xu[j,:] is the upper bound of box j 512 u, v, # old box bounds 513 u1,v1 # box in which the points are to be generated 514 ):

515 """ 516 If xnew contains points outside the box [u,v], the box is made larger 517 such that it just contains all new points; moreover, if necessary, 518 [u,v] is made larger to contain [u1,v1] 519 the box bounds xl and xu of the boxes on the boundary and the volume 520 measures small of these boxes are updated if necessary 521 522 Output: 523 ------- 524 xl updated version of xl 525 xu updated version of xu 526 u,v updated box bounds (the old box is contained in the new one) 527 528 And: 529 small updated version of small 530 """ 531 (nx,n) = xl.shape 532 uold = u.copy() 533 vold = v.copy() 534 535 u = numpy.minimum(numpy.min(xnew, axis=0), u, u1) 536 v = numpy.maximum(numpy.max(xnew, axis=0), v, v1) 537 538 i1 = find(u<uold) 539 i2 = find(v>vold) 540 541 ind = numpy.array([]) 542 for j in xrange( len(i1) ): 543 j1 = find(xl[:,i1[j]]==uold[i1[j]]) 544 ind = numpy.append(ind, j1) 545 xl[j1,i1[j]] = u[i1[j]] 546 547 for j in xrange( len(i2) ): 548 j2 = find(xu[:,i2[j]] == vold[i2[j]] ) 549 ind = numpy.append(ind, j2) 550 xu[j2,i2[j]] = v[ i2[j] ] 551 552 if len(i1) + len(i2): # At least one of the bounds was changed 553 self.small = self._calc_Smallness( xu, xl, v, u ) 554 555 return (xl, xu, u, v)

556 557

558 - def _selectIndexByVariance(self, x ):

559 """ 560 Choose the index i such that the variance of 561 x[:,i]/(self.v[i]-self.u[i]) is maxmal. 562 """ 563 variance = numpy.zeros(self.n) 564 for i in xrange(self.n): 565 variance[i] = std( x[:,i]/(self.v[i]-self.u[i]) ) 566 return numpy.argmax( variance )

567 568

569 - def _calc_Lambda(self, f1, f2):

570 """ 571 Calculate the lambda for splitting the box 572 """ 573 if f1 <= f2: 574 return self.goldenSection 575 else: 576 return 1.0 - self.goldenSection

577 578

579 - def _calc_SplitPoint(self, j, x1, x2):

580 """ 581 Compute the split point 582 """ 583 _lambda = self._calc_Lambda(self.func(self.x[j,:]), 584 self.func(self.x[j+1,:]) 585 ) 586 return _lambda*x1 + (1.0 - _lambda)*x2

587 588

589 - def _split(self, 590 x, # the rows are a set of points. 591 f, # f[j] contains the function value. 592 df, # its variation. 593 xl0, # vector of lower bounds of the box. 594 xu0, # vector of upper bounds of the box. 595 nspl # nspl[i] is the number of splits the box [xl0,xu0]. 596 # has already undergone along the ith coordinate. 597 ):

598 """ 599 Splits a box [xl0,xu0] contained in a bigger box [u,v] such that each 600 of the resulting boxes contains just one point of a given set of points 601 602 Output: 603 ------ 604 xl xl[j,:] is the lower bound of box j 605 xu xu[j,:] is the upper bound of box j 606 x x[j,:] is the point contained in box j 607 f f[j] contains the function value at x[j,:]. 608 df df[j] the uncertainty of f[j]. 609 nsplit nsplit[j,i] is the number of times box j has been split in the 610 i-th coordinate 611 small small[j] is an integer-valued logarithmic volume measure of box j 612 """ 613 (nx,n) = x.shape 614 if nx == 1: 615 small = self._calc_Smallness( xu0, xl0, self.v, self.u ) 616 return (xl, xu, x, f, df, nspl, small ) 617 618 elif nx == 2: 619 i = numpy.argmax( abs(x[0,:]-x[1,:]) / (self.v-self.u) ) 620 621 if self.constraint is None: 622 _lambda = self._calc_Lambda(f[0], f[1]) 623 ymid = _lambda*x[0,i] + (1.0 - _lambda)*x[1,i] 624 else: 625 ymid = 0.5*(x[0,i] + x[1,i]) 626 627 xl = numpy.zeros( (2,len(xl0) ) ) 628 xu = numpy.zeros( (2,len(xu0) ) ) 629 630 xl[0,:] = xl0; xl[1,:] = xl0 631 xu[0,:] = xu0; xu[1,:] = xu0 632 633 if x[0,i] < x[1,i]: 634 xu[0,i] = ymid; xl[1,i] = ymid 635 else: 636 xl[0,i] = ymid; xu[1,i] = ymid 637 638 nsplit = numpy.zeros( (2, len(nspl) ), 'int' ) 639 nsplit[0,:] = nspl 640 nsplit[0,i] += 1 641 nsplit[1,:] = nsplit[0,:] 642 small = self._calc_Smallness( xu, xl, self.v, self.u ) 643 644 return (xl, xu, x, f, df, nsplit, small ) 645 646 647 # Choose the index i such that the variance of 648 # x[:,i]/(self.v[i]-self.u[i]) is maxmal. 649 i = self._selectIndexByVariance( x ) 650 651 # Sort the points 652 y = self.rsort( x[:,i] ) 653 654 j = numpy.argmax( y[1:len(y)]-y[0:len(y)-1] ) 655 if self.constraint is None: 656 ymid = self._calc_SplitPoint(j, y[j], y[j+1]) 657 else: 658 ymid = 0.5*(y[j]+y[j+1]) 659 660 ind1 = find( x[:,i]<ymid ) 661 ind2 = find( x[:,i]>ymid ) 662 663 xl = numpy.zeros( (2,len(xl0)) ) 664 xu = numpy.zeros( (2,len(xu0)) ) 665 xl[0,:] = xl0; xl[1,:] = xl0; xl[1,i] = ymid 666 xu[0,:] = xu0; xu[1,:] = xu0; xu[0,i] = ymid 667 668 nsplit = numpy.zeros( (2,len(xl0)), 'int' ) 669 nsplit[0,:] = nspl 670 nsplit[0,i] += 1 671 nsplit[1,:] = nsplit[0,:] 672 673 npoint = numpy.array( [len(ind1), len(ind2)] ) 674 ind = -1*numpy.ones( (2, max(len(ind1),len(ind2))),'int' ) 675 ind[ 0,0:len(ind1) ] = ind1 676 ind[ 1,0:len(ind2) ] = ind2 677 678 nboxes = 1 679 [maxpoint,j] = max_(npoint) 680 while int(maxpoint) > 1: 681 682 ind0 = ind[j,find(ind[j,:]>=0)] 683 684 # Choose the index i such that the variance of 685 # x[ind0:,i]/(self.v[i]-self.u[i]) is maxmal. 686 i = self._selectIndexByVariance( x[ind0,:] ) 687 688 # Sort the points 689 y = self.rsort( x[ind0,i] ) 690 d = y[ 1:len(y) ] - y[ 0:len(y)-1 ] 691 if len(d) ==0: 692 break 693 694 k = numpy.argmax( d ) 695 if self.constraint is None: 696 ymid = self._calc_SplitPoint(k, y[k], y[k+1]) 697 else: 698 ymid = 0.5*(y[k] + y[k+1]) 699 700 ind1 = ind0[ find( x[ind0,i]<ymid ) ] 701 ind2 = ind0[ find( x[ind0,i]>ymid ) ] 702 703 nboxes += 1 704 705 xl = numpy.append(xl, [ xl[j,:] ], axis=0 ) 706 xu = numpy.append(xu, [ xu[j,:] ], axis=0 ) 707 xu[j,i] = ymid 708 xl[nboxes,i] = ymid 709 710 nsplit[j,i] += 1 711 nsplit = numpy.append(nsplit, [ nsplit[j,:] ], axis=0 ) 712 npoint[j] = len(ind1) 713 npoint = numpy.append(npoint, len(ind2) ) 714 715 ind[j,0:ind.shape[1] ] = -1 716 ind[j,0:len(ind1) ] = ind1 717 ind = numpy.append(ind,[-1*numpy.ones(ind.shape[1],'int')], axis=0) 718 ind[nboxes,0:len(ind2)] = ind2 719 720 [maxpoint,j] = max_(npoint) 721 722 # ---------------------------------- 723 x = x[ ind[:,0], :] 724 f = f[ ind[:,0] ] 725 df = df[ ind[:,0] ] 726 small = self._calc_Smallness( xu, xl, self.v, self.u ) 727 728 return ( xl, xu, x, f, df, nsplit, small )

729 730

731 - def _notnanMaxMin(self, notnan):

732 """ This function should be merged latter """ 733 if len(notnan) !=0 : 734 fmin = min( self.f[notnan] ) 735 fmax = max( self.f[notnan] ) 736 else: 737 fmin = 1 738 fmax = 0 739 return (fmax,fmin)

740 741

742 - def _nanMaxMin(self):

743 """ 744 Find the max and min function values, but ignore nan elements 745 """ 746 fmin = numpy.nanmin(self.f) 747 fmax = numpy.nanmax(self.f) 748 return (fmax,fmin)

749 750

751 - def _filterInput(self, 752 x, # the rows of x are a set of points 753 f, # f[j,:] is function value of x[j,:] & its uncertainty 754 df 755 ):

756 """ 757 Checks whether there are any duplicates among the points given by the 758 rows of x, throws away the duplicates and computes their average 759 function values and an estimated uncertainty 760 761 See details in Algorithm Step 2 762 763 Output: 764 ------ 765 x updated version of x (possibly some points have been deleted) 766 f updated version of f ( the average function value ). 767 df updated version of df (the estimated uncertainty pertaining to x ) 768 nxFval nxFval[j] is the number of times the row x[j,:] appeared 769 in the input version of x 770 t t[j] is nxFval[j] times the variance of the function values 771 measured for point x[j,:] 772 """ 773 (nx,n) = x.shape 774 # define np and t 775 np = ivector(nx) 776 t = ivector(nx) 777 778 i = 0 779 while i < nx: 780 j = i+1 781 ind = [] 782 while j < nx: 783 if sum( x[i,:]-x[j,:] ) == 0: 784 ind.append(j) 785 j += 1 786 787 if len(ind) !=0 : 788 ind = ind.insert(0,i) 789 ind1 = find( numpy.isnan( f[ind,0] ) ) 790 791 if len(ind1) < len(ind): 792 ind[ind1] = [] 793 np[i] = len( ind ) 794 fbar = sum( f[ind])/np[i] 795 t[i] = sum( (f[ind]-fbar)**2 ) 796 f[i] = fbar 797 df[i] = numpy.sqrt( (sum( df[i]**2 ) + t[i])/np[i] ) 798 else: 799 np[i] = 1 800 t[i] = 0 801 802 # More test here 803 x[ind[1:end],:] = [] 804 f[ind[1:end]] = [] 805 df[ind[1:end]] = [] 806 nx = x.shape[0] 807 else: 808 np[i] = 1 809 t[i] = 0 810 811 i += 1 812 813 return [x,f,df,np,t]

814 815

816 - def _calcSnnByInd(self, i):

817 """ Computes safeguarded nearest neighbors 818 819 Computes a safeguarded set of nearest neighbors to a point 820 x0 = self.x[i,:] for each coordinate, the nearest point differing 821 from x0 in the i-th coordinate by at least 1 is chosen 822 self.x are the points from which the neighbors are to 823 be chosen (possibly x0 is among these points) 824 i: The index of the point around which the fit is to be computed 825 826 Output: 827 ------- 828 near Vector pointing to the nearest neighbors (i.e. to the rows of x) 829 d Maximal distance between x0 and a point of the set of 830 nearest neighbors 831 """ 832 x0 = self.x[i,:] # Point for which the neighbors are to be found 833 d = numpy.sqrt( numpy.sum( (self.x- 834 numpy.array([x0]).repeat(self.nx,axis=0) 835 )**2, 836 axis=1 837 ) 838 ) 839 [d1,ind] = sort(d) 840 841 # Eliminate x0 if it is in the set 842 if not d1[0]: 843 ind = numpy.delete(ind, 0) 844 845 near = ivector() 846 for i in xrange(self.n): 847 j = min(find(abs(self.x[ind,i]-x0[i])-self.dx[i] >= 0)) 848 near = numpy.append(near, ind[j]) 849 ind = numpy.delete(ind, j) 850 851 j = 0 852 while near.size < self.snn: # snn: The number of neighbors to be found 853 if len(find(near==ind[j]))==0 and \ 854 max( abs(x0-self.x[ind[j],:]) - self.dx ) >= 0: 855 near = numpy.append( near, ind[j] ) 856 j +=1 857 858 [d,ind] = sort( d[near] ) 859 860 return ( near[ind], max(d) )

861 862

863 - def _updateFromNewPoints(self, 864 xl, # lower bound of the old boxes 865 # its variation and other parameters 866 xu, # upper bounds of the old boxes 867 # its variation and other parameters 868 xnew, # rows contain new points 869 fnew, # new function values and 870 dfnew # their variations 871 ):

872 """ 873 Updates the box parameters when a set of new points and their function 874 values are added, i.e., the boxes containing more than one point are 875 split and the nearest neighbors are computed or updated. 876 877 Note: 878 ----- 879 self.u1, 880 self.v1: box in which the new points are to be generated 881 self.x: rows contain the old points 882 self.f: f[j] contains the function value at x[j,:] 883 884 Output: 885 ------- 886 xl,xu updated version of xl,xu (including new boxes) 887 inew pointer pointing to the new boxes and boxes whose 888 nearest neighbors have changed 889 890 And possibly updated version of: 891 self.u, self.v( updated box bounds such that all new points are in box) 892 self.fnan updated version of fnan (if a function value 893 was found for a point in the new iteration) 894 self.near, self.small, self.d, self.nsplit 895 self.nxFval, self.t, self.x, self.f 896 """ 897 # The number of points from the previous iteration 898 nxold = self.nx 899 nxnew = xnew.shape[0] 900 inew = ivector() 901 902 # Do we need Case nxold==0??? 903 # If any of the new points are already among old points, they are 904 # thrown away and function value and its uncertainty are updated 905 _del = ivector() 906 for j in xrange(nxnew): 907 908 i=self._sumFind(abs(numpy.array([xnew[j,:]]).repeat(nxold,axis=0)- 909 self.x), 910 0 911 ) 912 if len(i) != 0: 913 914 if len( find(self.fnan==i) )==0 and \ 915 numpy.isfinite( self.f[i] ) : 916 917 if not numpy.isnan( fnew[j] ): 918 self.nxFval[i] += 1 919 delta = fnew[j] - self.f[i] 920 self.f[i] += delta/self.nxFval[i] 921 self.t[i] += delta*(fnew[j]-self.f[i]) 922 self.df[i] = numpy.sqrt( self.df[i]**2 + \ 923 ( delta*(fnew[j]-self.f[i])+ \ 924 fnew[j]**2-self.df[i]**2)/self.nxFval[i] 925 ) 926 inew = numpy.append(inew,i) 927 928 else: # point i had NaN function value 929 930 if not numpy.isnan( fnew[j] ): 931 self.f[i] = fnew[j] 932 self.df[i] = dfnew[j] 933 inew = numpy.append(inew,i) 934 if len(self.fnan) != 0: 935 ii = find(self.fnan==i) 936 self.fnan = removeByInd(self.fnan, ii) 937 938 _del = numpy.append( _del, j ) 939 940 941 if len(_del) != 0 : 942 xnew = numpy.delete( xnew, _del, axis=0 ) 943 fnew = numpy.delete( fnew, _del, axis=0 ) 944 dfnew = numpy.delete(dfnew, _del, axis=0 ) 945 946 nxnew = xnew.shape[0] 947 if nxnew==0: 948 inew = sort(inew) 949 return (xl, xu, inew) 950 951 # Filter the new recommended evaluation points 952 [xnew, fnew0, dfnew0, npnew, tnew] = \ 953 self._filterInput(xnew, fnew, dfnew ) 954 nxnew = xnew.shape[0] 955 956 # update current number of points 957 self.nx = nxold + nxnew 958 959 # update the new bounds if necessary 960 if numpy.sum(numpy.minimum(numpy.min(xnew,axis=0),self.u)<self.u) or \ 961 numpy.sum(numpy.maximum(numpy.max(xnew,axis=0),self.v)>self.v) or \ 962 numpy.sum( numpy.minimum(self.u,self.u1) < self.u ) or \ 963 numpy.sum( numpy.maximum(self.v,self.v1) > self.v ): 964 [xl, xu, self.u, self.v] = \ 965 self._newBounds(xnew, xl, xu, self.u, self.v, self.u1, self.v1) 966 967 self.x = numpy.append(self.x, xnew, axis=0) 968 inew = numpy.append(inew, range(nxold,self.nx) ) 969 970 for i in range(nxold,self.nx): 971 self.f = numpy.append(self.f, fnew0[i-nxold] ) 972 self.df = numpy.append(self.df, dfnew0[i-nxold] ) 973 974 self.nxFval = numpy.append(self.nxFval, npnew) 975 self.t = numpy.append(self.t, tnew) 976 977 par = ivector( nxnew ) 978 for j in xrange(nxnew): 979 ind=self._sumFind(within(xl, 980 numpy.array([xnew[j,:]]).repeat(nxold, 981 axis=0), 982 xu), 983 self.n 984 ) 985 if len(ind) !=0 : 986 par[j] = ind[ numpy.argmin(self.small[ind]) ] 987 988 par1 = self.rsort(par) 989 inew = numpy.append(inew, par1 ) 990 991 # Add more space for xl and xu. Check the size is nxnew ??? 992 xl = numpy.append(xl, 993 numpy.array([vector(self.n)]).repeat(nxnew, axis=0), 994 axis=0 995 ) 996 xu = numpy.append(xu, 997 numpy.array([vector(self.n)]).repeat(nxnew, axis=0), 998 axis=0 999 ) 1000 self.nsplit=numpy.append(self.nsplit, 1001 numpy.array([ivector(self.n)]).repeat(nxnew, 1002 axis=0), 1003 axis=0 1004 ) 1005 self.small =numpy.append(self.small, numpy.zeros(nxnew) ) 1006 for l in xrange( len(par1) ): 1007 j = par1[l] 1008 ind = find(par==j) + nxold 1009 idx = numpy.append( numpy.array([j]), ind, axis=0 ) 1010 [xl0, xu0, x0, f0, df0, nsplit0, small0] = \ 1011 self._split(self.x[idx,:], 1012 self.f[idx], 1013 self.df[idx], 1014 xl[j,:], 1015 xu[j,:], 1016 self.nsplit[j,:] 1017 ) 1018 k = self._sumFind( x0 == \ 1019 numpy.array([self.x[j,:]]).repeat( 1020 x0.shape[0],axis=0), 1021 self.n 1022 ) 1023 if len(k) > 1: 1024 k = k[0] # Choose the first one 1025 1026 xl[j,:] = xl0[k,:] 1027 xu[j,:] = xu0[k,:] 1028 self.nsplit[j,:] = nsplit0[k,:] 1029 self.small[j] = small0[ int(k) ] 1030 1031 for k in xrange(len(ind)): # number of pts in box [xl[j,:],xu[j,:]] 1032 k1 = ind[k] 1033 k2 = self._sumFind(x0 == \ 1034 numpy.array([self.x[k1,:]]).repeat( 1035 x0.shape[0],axis=0), 1036 self.n 1037 ) 1038 if len(k2) > 1: 1039 k2 = k2[ numpy.argmin(self.small[k2]) ] 1040 1041 # More test Later 1042 xl[k1,:] = xl0[k2,:] 1043 xu[k1,:] = xu0[k2,:] 1044 self.nsplit[k1,:] = nsplit0[k2,:] 1045 self.small[k1] = small0[k2] 1046 1047 # Find the current max and min function values. 1048 (fmax, fmin ) = self._nanMaxMin() 1049 1050 if self.nx >= self.snn+1 and fmin < fmax: 1051 for j in range(nxold,self.nx): 1052 [nn,dd] = self._calcSnnByInd( j ) 1053 self.d = numpy.append(self.d, float(dd) ) 1054 self.near = numpy.append(self.near, [nn], axis=0 ) 1055 1056 for j in xrange(nxold): 1057 xx = ( numpy.array([self.x[j,]]).repeat(nxnew,axis=0)-xnew )**2 1058 if min( numpy.sqrt( numpy.sum( xx, axis=1 ) ) ) < self.d[j]: 1059 [ self.near[j,:], self.d[j] ] = self._calcSnnByInd( j ) 1060 inew = numpy.append(inew, j) 1061 1062 inew = self.rsort(inew) 1063 else: 1064 self.near = ivector() 1065 self.d = inf*ones(1,self.nx) 1066 1067 return (xl, xu, inew)

1068 1069 1070 1071 # -----------------------------------------------------------------------

1072 - def _replace_nan(self, 1073 inew # vector pointing to the new boxes and boxes 1074 # whose nearest neighbors have changed 1075 ):

1076 """ 1077 Replaces the function values NaN of a set of points 1078 by a value determined by their nearest neighbors with finite function 1079 values, with a safeguard for the case that all neighbors 1080 have function value NaN 1081 """ 1082 notnan = range(1, self.f.size) 1083 if len(self.fnan): 1084 notnan[self.fnan] = [] 1085 1086 [fmax,imax] = max_( self.f ) 1087 fmin = min(self.f) 1088 1089 dfmax = self.df[imax] 1090 1091 for j in range( len(self.fnan) ): 1092 1093 l = self.fnan[j] 1094 if len(find(inew==l)) != 0: 1095 # A substitute function value is only computed for new points 1096 # and for points whose function values have changed 1097 ind = self.near[l,:] 1098 1099 # Eliminate neighbors with function value NaN 1100 for i in xrange( len(ind) ): 1101 if len( find(self.fnan==ind[i]) ) != 0: 1102 ind = numpy.delete( ind, numpy.s_[ i ], axis=0 ) 1103 1104 # Updated self.f 1105 if len(ind)==0: 1106 self.f[l] = fmax + 1.e-3*(fmax-fmin) 1107 self.df[l] = dfmax 1108 else: 1109 [fmax1,k] = numpy.max( self.f[ind] ) 1110 fmin1 = numpy.min( self.f[ind] ) 1111 self.f[l] = fmax1 + 1.e-3*(fmax1-fmin1) 1112 self.df[l] = self.df[ ind[k] ]

1113 1114 1115 1116 # ------------------------------------------------------------------------

1117 - def _genClass5Points(self, 1118 x, # The points chosen by Snobfit (can be empty) 1119 nreq # The number of points to be generated 1120 ):

1121 """ 1122 Generates nreq recommended evalution points of Class 5 1123 in [self.u1, self.v1] 1124 1125 [self.u,self.v]: Bounds of the box in which the points to be generated 1126 1127 Output: 1128 ------- 1129 x1 The rows are the requested points 1130 x1 is of dimension nreq x n (n = dimension of the problem ) 1131 """ 1132 nx1 = 100*nreq 1133 1134 xnew = numpy.round( numpy.array([self.u1]).repeat(nx1, axis=0) + \ 1135 rand(nx1,self.n)* \ 1136 numpy.array([self.v1-self.u1]).repeat(nx1, axis=0) 1137 ) 1138 nx = x.shape[0] 1139 if nx: 1140 d = vector(nx1) 1141 for j in xrange(nx1): 1142 xnew[j,:] = self.round(xnew[j,:], self.u1, self.v1) 1143 d[j]=numpy.min(numpy.sum((x - numpy.array([xnew[j,:]] 1144 ).repeat(nx, axis=0) )**2, 1145 axis=0) 1146 ) 1147 ind = find( d==0 ) 1148 1149 if len(ind) != 0: 1150 xnew[ind,:] = [] # Fix later 1151 d[ind] = [] 1152 1153 x1 = numpy.array( [] ) 1154 nx1 = xnew.shape[0] 1155 1156 else: # TEST Later 1157 x1 = xnew[0,:] 1158 xnew = numpy.delete( xnew, numpy.s_[ 0 ], axis=0 ) 1159 nx1 -= 1 1160 d = numpy.sum( (xnew - numpy.array([x1]).repeat(nx1, axis=0) )**2, 1161 axis=0 1162 ) 1163 nreq -= 1 1164 1165 for j in xrange(nreq): 1166 i = numpy.argmax(d) 1167 y = xnew[i,:] 1168 if j == 0 and len(x1)==0: 1169 x1 = numpy.mat([y]) 1170 else: 1171 x1 = numpy.append(x1, [y], axis=0) 1172 1173 xnew = numpy.delete( xnew, numpy.s_[ i ], axis=0 ) 1174 d = numpy.delete( d, numpy.s_[ i ], axis=0 ) 1175 nx1 -= 1 1176 1177 d1 = numpy.sum( ( xnew-numpy.array([y]).repeat(nx1, axis=0) )**2, 1178 axis=1 1179 ) 1180 d = numpy.minimum( d,d1 ) 1181 1182 return numpy.asarray(x1)

1183 1184 1185 #------------------------------------ 1186 # Step 1 1187 #------------------------------------

1188 - def update_uv(self):

1189 """ Update the search bounds [self.u, self.v] 1190 1191 The search box [self.u, self.v] is updated if the user changed it. 1192 The vectors self.u and self.v are defined such that [self.u, self.v] 1193 is the smallest box containing [u', v'], all new input points and, 1194 in the case of a continuation call to Snobfit, also the box [u1, v1] 1195 from the previous iteration. 1196 [self.u, self.v] is considered to be the box to be explored and 1197 a box tree of [self.u, self.v] is generated; 1198 1199 Note: All suggested new evaluation points are in [u', v']. 1200 """ 1201 if len(self.x) != 0: 1202 self.u = numpy.minimum( numpy.min(self.x,axis=0), self.u1 ) 1203 self.v = numpy.maximum( numpy.max(self.x,axis=0), self.v1 )

1204 1205 1206 #------------------------------------ 1207 # Step 2 1208 #------------------------------------

1209 - def update_input(self):

1210 """ 1211 Throw out duplicates among the points and 1212 compute the mean function value and deviation ( ie., self.f, self.df ) 1213 """ 1214 (self.x, self.f, self.df, self.nxFval, self.t) = \ 1215 self._filterInput(self.x, self.f, self.df)

1216 1217 1218 #------------------------------------ 1219 # Step 3 1220 #------------------------------------

1221 - def branch(self):

1222 """ 1223 All current boxes containing more than one point are split 1224 according to the algorithm described in Section 5. 1225 The smallness (also defined in Section 5) is computed for these boxes. 1226 If [self.u, self.v] is larger than [uold, vold] in a continuation call, 1227 the box bounds and the smallness are updated for the boxes 1228 for which this is necessary. 1229 """ 1230 if self.nx: 1231 [self.xl,self.xu,self.x,self.f,self.df,self.nsplit,self.small] = \ 1232 self._split(self.x, 1233 self.f, 1234 self.df, 1235 self.u, 1236 self.v, 1237 numpy.zeros( self.n ) 1238 ) 1239 else: 1240 self.xl = vector() 1241 self.xu = vector() 1242 self.nsplit = ivector() 1243 self.small = ivector()

1244 1245 1246 #------------------------------------ 1247 # Step 4 1248 #------------------------------------

1249 - def calc_snn(self):

1250 """ Compute the safeguarded nearest neighbors. 1251 1252 If we have |X| < n + dn + 1 for the current set X of evaluation points, 1253 or if all function values != NaN (if any) are identical, go to Step 11. 1254 Otherwise, for each new point x a vector pointing to n + dn 1255 safeguarded nearest neighbors is computed. The neighbor lists of 1256 some old points are updated in the same way if necessary. 1257 1258 Note: we denote the X the set of points for which the objective 1259 function has already been evaluated at some stage of SNOBFIT. 1260 """ 1261 for i in xrange( self.nx ): 1262 [ self.near[i,:], self.d[i] ] = self._calcSnnByInd(i)

1263 1264 1265 #------------------------------------ 1266 # Step 5 1267 #------------------------------------

1268 - def calc_fdf_nan(self):

1269 """ Calculate the fictitious values f and df for nan points 1270 1271 For any new input point x with function value NaN (which means that 1272 the function value could not be determined at that point) and for all 1273 old points marked infeasible whose nearest neighbors have changed, 1274 we define fictitious values f and df as specified in Section 3. 1275 1276 Later we test an objective function with nan values 1277 """ 1278 self.fnan = find( isnan( self.f ) ) 1279 if not isEmpty(self.fnan): 1280 self._replace_nan( range(0,self.nx) )

1281 1282

1283 - def _fitLocalModel(self, j ):

1284 """ 1285 Computes a local fit around the point x0 = x[j,:] 1286 and minimizes it on a trust region 1287 j: The index of the point around which the fit is to be computed 1288 1289 See details in Algorithm Step 6 1290 1291 Output: 1292 ------- 1293 y Estimated minimizer in the trust region 1294 fy Its estimated function value 1295 g Estimated gradient for the fit 1296 sigma sigma = norm(A*g-b)/sqrt(K-n), 1297 where A and b are the coefficients resp. right hand side of the 1298 fit, n is the dimension and K the number of nearest neighbors 1299 considered(estimated standard deviation of the model errors) 1300 """ 1301 (x0, f0, df0) = ( self.x[j,:], self.f[j], self.df[j] ) 1302 D = self._calc_DMat(df0) 1303 x1 = self.x[ self.near[j,:],:] 1304 K = x1.shape[0] 1305 S = x1 - numpy.array( [x0] ).repeat(K, axis=0) 1306 1307 d = numpy.maximum(0.5*numpy.max(abs(S), axis=0), self.dx) 1308 1309 # Set the Q_k = (x_k-x) * D * (x _k -x ) 1310 Q = vector(K) 1311 for k in xrange(K): 1312 Q[k] = dot3( S[k,:],D,toCol(S[k,:]) ) + self.df[ self.near[j,k] ] 1313 1314 # Now the A_ki = (x_i-x_i^k)/Q_k 1315 A = S / numpy.array( toCol(Q) ).repeat(self.n, axis=1) 1316 1317 # Now the b_k = (f-f_k) /Q_k 1318 b = toCol ( ( self.f[ self.near[j,:] ] - f0) / Q ) 1319 1320 # Reduced SVD 1321 [U, Sigma, V] = svd(A) 1322 1323 # To enforce a condition number < 10^-4 1324 Sigma = numpy.maximum(Sigma, 1.e-4*Sigma[0] ) 1325 1326 # To calc the Estimated gradient 1327 g = dot3( numpy.dot(V, numpy.diag( 1.0 / Sigma )), U.T, b) 1328 1329 # Calc the estimated standard deviation of the model errors 1330 residual = numpy.array( (A*toCol(g) - b) ) **2 1331 sigma = numpy.sqrt( numpy.sum( residual )/(K-self.n) ) 1332 1333 # In the trust region, calc the Estimated minimizer 1334 pl = numpy.maximum(-d, self.u-x0) # lower bounds of p 1335 pu = numpy.minimum( d, self.v-x0) # upper bounds of p 1336 p = vector(self.n) # p = y-x0 1337 for i in xrange(self.n): 1338 p[i] = self.quadMin( float(sigma*D[i,i]), g[i], pl[i], pu[i] ) 1339 1340 y = self.round(x0+p, self.u, self.v) 1341 1342 nc = 0 1343 while min(numpy.max(abs(self.x - \ 1344 numpy.array([y]).repeat(self.nx, axis=0)) - \ 1345 numpy.array([self.dx]).repeat(self.nx, axis=0), 1346 axis=1))<0 and nc < 5: 1347 p = pl + (pu-pl)*rand(self.n) 1348 y = self.round(x0+p, self.u, self.v) 1349 nc += 1 1350 1351 # Calc the Estimated Function Value 1352 fy = self._estimate(self.f[j], self.df[j], y, x0, g, sigma) 1353 1354 # Return g as row vector 1355 g = toRow(g) 1356 1357 return (y, fy, g, sigma)

1358 1359 1360 #------------------------------------ 1361 # Step 6 1362 #------------------------------------

1363 - def localFit(self):

1364 """ Create local surrogate models. 1365 1366 Local surrogate models are created by linear least squares fit. 1367 Local fits (as described in Section 3) around all new points and 1368 all old points with changed nearest neighbors are computed and 1369 the potential points y of Class 2 and 3 and their estimated 1370 function values fy are determined as in Section 4. 1371 """ 1372 for j in xrange(self.nx): 1373 self.y[j,:], self.fy[j], self.gradient[j,:], self.sigma[j] = \ 1374 self._fitLocalModel( j )

1375 1376

1377 - def _nanRequest(self, a):

1378 b = numpy.array([nan,5]).repeat(self.nreq, axis=0) 1379 # ( [nan, 5], self.nreq ) 1380 return numpy.append(a,b, axis=1)

1381 1382

1383 - def returnNanRequest(self):

1384 # We don't test this part carefully 1385 # Test this part latter. Fix me latter 1386 x5 = self._genClass5Points(self.x, self.nreq) 1387 self.request = self._nanRequest(x5) 1388 if len(self.x) != 0 : 1389 [ fbest, ind] = min_( self.f ) 1390 xbest = self.x[ind,:] 1391 1392 else: 1393 xbest = nan*numpy.ones(self.n) 1394 fbest = inf 1395 1396 if self._isNotEnoughReq(): 1397 self._warning() 1398 1399 return (self.request, xbest, fbest )

1400 1401

1402 - def _extend(self, a, n):

1403 """ 1404 Extend the Array a about n element(s). 1405 """ 1406 if len(a.shape)==1: nadd = numpy.zeros( n ) 1407 else: nadd = numpy.zeros( (n, a.shape[1]) ) 1408 return numpy.append(a, numpy.array(nadd), axis=0)

1409 1410

1411 - def updateWorkspace(self):

1412 """ update the current work space """ 1413 oldxbest = self.xbest 1414 1415 [self.xl, self.xu, self.inew] = \ 1416 self._updateFromNewPoints(self.xl, 1417 self.xu, 1418 self.xnew, 1419 self.fnew, 1420 self.dfnew 1421 ) 1422 1423 if not isEmpty(self.near): 1424 ind = find( isnan(self.f) ) 1425 self.fnan = numpy.append(self.fnan, ind) 1426 if not isEmpty(self.fnan): 1427 self._replace_nan( self.inew ) 1428 1429 # Extend more space for y, g, and sigma. 1430 nadd = max(self.inew) - self.y.shape[0] + 1 1431 if nadd > 0: 1432 self.y = self._extend( self.y, nadd ) 1433 self.fy = self._extend( self.fy, nadd ) 1434 self.sigma = self._extend( self.sigma, nadd ) 1435 self.gradient = self._extend( self.gradient, nadd ) 1436 1437 for j in self.inew: 1438 self.y[j], self.fy[j], self.gradient[j,:], self.sigma[j] = \ 1439 self._fitLocalModel(j)

1440 1441 1442 1443 # -----------------------------------------------------------------

1444 - def _quadFit(self, j ):

1445 """ 1446 A quadratic model around the best point is fitted and minimized with 1447 minq over a trust region 1448 j is the index of the best point 1449 1450 Output: 1451 ------- 1452 y minimizer of the quadratic model around the best point 1453 fy its estimated function value 1454 """ 1455 (x0, f0) = (self.x[j,:], self.f[j]) 1456 1457 K = min( self.nx-1, self.n*(self.n+3) ) 1458 # K = min( self.nx-1, self.n*(self.n+3)/2 + self.dn ) 1459 distance=numpy.sum((self.x-numpy.array([x0]).repeat(self.nx,axis=0))**2, 1460 axis=1 1461 ) 1462 ind = numpy.argsort(distance)[1:K+1] 1463 1464 d = numpy.max(abs( self.x[self.near[j,:],:] - \ 1465 numpy.array([x0]).repeat(self.snn, axis=0) ), 1466 axis=0 1467 ) 1468 d = numpy.maximum(d, self.dx) 1469 1470 # S^k = x^k - x0 ( xbest ) 1471 S = self.x[ind,:] - numpy.array( [x0] ).repeat(K, axis=0) 1472 1473 R = triu( qr(S) ) 1474 L = inv(R).T 1475 1476 # Scaling Matrix H which esures affine invariance of fitting procedure 1477 # Note: 1) S^T * H * S = || R^-T * S ||^2 1478 # 2) The exponent 3/2 in the error term reflects the expected 1479 # cubic approx. order of the quadratic model 1480 sc = numpy.sum( numpy.dot(S,L.T)**2, axis=1) ** (3.0/2.0 ) 1481 b = (self.f[ind]-f0)/sc 1482 1483 xx = self.x[ind,:] - numpy.array( [x0] ).repeat(K, axis=0) 1484 xx2 = 0.5*xx*xx 1485 A = numpy.bmat( 'xx xx2') 1486 1487 for i in xrange(self.n-1): 1488 xx =( self.x[ind,i] - x0[i]) 1489 B = numpy.array( toCol(xx) ).repeat( self.n-i-1, axis=1) 1490 xx = B*(self.x[ind,i+1:self.n] - \ 1491 numpy.array( [x0[i+1:self.n]] ).repeat(K, axis=0) ) 1492 A = numpy.bmat('A xx') 1493 1494 xx = numpy.array( toCol(sc) ).repeat(A.shape[1], axis=1) 1495 A = A/xx 1496 1497 # In q, it includes the components of g and the upper triangle of G 1498 q = mldiv(A,toCol(b)) 1499 1500 # G is a symmetric matrix 1501 G = numpy.diag( q[self.n:2*self.n] ) 1502 k = 2*self.n 1503 for i in xrange( self.n-1 ): 1504 for j in xrange(self.n-i-1): 1505 G[i, j+i+1] = q[k] 1506 G[j+i+1, i] = q[k] 1507 k += 1 1508 1509 # The estimate gradient 1510 g = q[0:self.n] 1511 1512 # Solve the minimizer of the quadratic model around the best point x0 1513 if x0.shape[0] == 1: 1514 y = float( self.quadMin( float(G), 1515 float(g-numpy.dot(G,x0)), 1516 float(numpy.maximum(x0-d, self.u1)), 1517 float(numpy.minimum(x0+d, self.v1)) 1518 ) ) 1519 else: 1520 [y,info] = minq(g-numpy.dot(G, x0), 1521 G, 1522 numpy.maximum(x0-d, self.u1), 1523 numpy.minimum(x0+d, self.v1), 1524 ) 1525 y = self.round(y, self.u1, self.v1) 1526 1527 nc = 0 1528 while nc <10 and \ 1529 min(numpy.max( abs( self.x - \ 1530 numpy.array([y]).repeat(self.nx,axis=0) ) - \ 1531 numpy.array( [self.dx] ).repeat(self.nx, axis=0), 1532 axis=1))<-eps: 1533 1534 u1 = numpy.maximum(x0-d, self.u1) 1535 v1 = numpy.minimum(x0+d, self.v1) 1536 1537 y = u1 + numpy.dot( rand(1,self.n), (v1-u1) ) 1538 y = self.round(y, self.u1, self.v1) 1539 nc += 1 1540 1541 fy = self._quadEstimate(f0, y, x0, g, G) 1542 1543 return (y, fy )

1544 1545 1546 #------------------------------------ 1547 # Step 7 1548 #------------------------------------

1549 - def locQuadFit(self):

1550 """ 1551 The current best point xbest and the current best function value fbest 1552 in [self.u1, self.v1] are determined. If the objective function 1553 has not been evaluated in [self.u1, self.v1] yet, (i.e., n1 = 0), 1554 it mean that we doesn't generate the recommended evalution points of 1555 Class 1. And go to Step 8. 1556 1557 A local quadratic fit around self.xbest is computed as described 1558 in Section 3 and the point z of Class 1 is generated as described 1559 in Section 4. If such a point z was generated, let z be contained in 1560 the subbox [\underline{x}^j, \bar{x}^j ] (in the case that z belongs 1561 to more than one box, a box with minimal smallness is selected). 1562 if 1563 min_i( (\bar{x}^j_i - \underline{x}^j_i )/(v_i-u_i) ) > 1564 0.05 * max_i( (\bar{x}^j_i - \underline{x}^j_i )/(v_i-u_i) ), 1565 1566 the point z is put into the list of recommended evaluation points, 1567 and otherwise (if the smallest side length of the box 1568 [ \bar{x}_j , \underline{x}_j ] relative to [self.u, self.v] is too 1569 small compared to the largest one) we set self.J4 = {j}, i.e., 1570 a point of Class 4 is to be generated in this box later, 1571 which seems to be more appropriate for such a long, narrow box. 1572 This gives n1 recommended evaluation points (n1 = 0 or 1) of Class 1. 1573 """ 1574 if numpy.sum( within(self.u1,self.xbest,self.v1), axis=0) == self.n: 1575 j_Class1 = self.jbest 1576 else: 1577 j_Class1 = self.ind[ numpy.argmin(self.f[self.ind]) ] 1578 1579 [z, fz] = self._quadFit( j_Class1 ) 1580 z = self.round(z, self.u1, self.v1) 1581 1582 j = self._sumFind( within(self.xl, 1583 numpy.array([z]).repeat(self.nx, axis=0), 1584 self.xu), 1585 self.n 1586 ) 1587 if len(j) > 1: 1588 j = j[ numpy.argmin(self.small[j]) ] 1589 1590 if min( numpy.max(abs(self.x- \ 1591 numpy.array( [z] ).repeat(self.nx, axis=0) ) - \ 1592 numpy.array( [self.dx] ).repeat(self.nx, axis=0), 1593 axis=1) ) >= -eps: 1594 1595 if self._isMarkedForClass4(j,self.xu,self.xl,self.v,self.u): 1596 self.J4 = numpy.append(self.J4, j) 1597 else: 1598 self.request = numpy.array( [ self.newRequest(z,fz,1) ] )

1599 1600 1601 #------------------------------------ 1602 # Step 8 1603 #------------------------------------

1604 - def assignN4EachClass(self):

1605 """ Use random number to assign new evaluation points for each Class. 1606 1607 To generate the remaining self.globloc := nreq - n1 recommended 1608 evaluation points, let glob1 := floor(self.p*self.globloc) 1609 glob2 := ceil*(self.p*self.globloc). 1610 1611 Then a random number generator sets self.glob = glob1 with 1612 probability self.globloc*self.p - glob1 and self.glob=glob2 otherwise. 1613 Then self.globloc - self.glob points of Classes 2 and 3 together 1614 and self.glob points of Class 4 are to be generated. 1615 """ 1616 # the remaining recommended evaluation points 1617 self.globloc = self.nreq - self.request.shape[0] 1618 1619 glob1 = self.globloc * self.p 1620 glob2 = int( numpy.floor(glob1) ) 1621 if rand() < glob1 - glob2: self.glob = glob2 + 1 1622 else: self.glob = glob2

1623 1624

1625 - def _genPoints(self, k, local, classID = 2):

1626 """ a Helper Function """ 1627 j = 0 1628 sreq = self.request.shape[0] 1629 while sreq < self.nreq-self.glob and j < len(local): 1630 1631 i = local[ k[j] ] 1632 y = self.round( self.y[i,:], self.u1, self.v1) 1633 1634 l = self._sumFind( within(self.xl, 1635 numpy.array([y]).repeat(self.nx,axis=0), 1636 self.xu), 1637 self.n 1638 ) 1639 1640 if len(l) == 0: # To solve the empty bug 1641 j += 1 1642 continue 1643 1644 if len(l) > 1: 1645 l = l[ numpy.argmin(self.small[l]) ] 1646 1647 if self._isMarkedForClass4(l, self.xu,self.xl, self.v, self.u): 1648 self.J4 = numpy.append(self.J4, l) 1649 j += 1 1650 continue 1651 1652 yy = self.request[:,0:self.n] - \ 1653 numpy.array([y]).repeat(sreq, axis=0) 1654 dxy = numpy.array( [numpy.maximum(self.dy,self.dx)] 1655 ).repeat(sreq, axis = 0) 1656 if numpy.max( abs(y-self.x[l,:])-self.dx) >= -eps and \ 1657 ( not sreq or \ 1658 min( numpy.max(abs(yy)-dxy, axis=1) ) >= -eps): 1659 1660 if sum(y==self.y[i,:]) < self.n: 1661 f = self._estimate(self.f[i], 1662 self.df[i], 1663 y, 1664 self.x[i,:], 1665 self.gradient[i,:], 1666 self.sigma[i] 1667 ) 1668 else: 1669 f = self.fy[i] 1670 1671 self._addRequest( [self.newRequest(y, f, classID)] ) 1672 1673 sreq = self.request.shape[0] 1674 j += 1

1675 1676 1677 1678 # ---------------------------------------------------------------------

1679 - def _calcLocalPoints(self, 1680 ind=None #pointer to the boxes to be considered 1681 ):

1682 """ 1683 Computes a pointer to all `local' points (i.e. points whose neighbors 1684 have `significantly larger' function values) 1685 1686 Output: 1687 ------- 1688 local vector containing the indices of all local points 1689 nonLocal vector containing the indices of all nonlocal points 1690 """ 1691 if ind==None: 1692 ind = range(0,self.nx) 1693 1694 jj = ivector() 1695 local = ivector() 1696 nonLocal = ind 1697 1698 for j in xrange( len(ind) ): 1699 ( fmax,fmin ) = self._maxmin( self.f[ self.near[ind[j],:] ] ) 1700 if self.f[ ind[j] ] < fmin-0.2*(fmax-fmin): 1701 local = numpy.append(local, ind[j] ) 1702 jj = numpy.append(jj, j) 1703 1704 if len(jj) != 0: 1705 nonLocal = removeByInd( nonLocal, jj ) 1706 1707 return (local,nonLocal)

1708 1709

1710 - def _constructPoints(self, i):

1711 """ 1712 For a box [self.xl, self.xu] containing a point x, a point y in the 1713 intersection of [self.xl, self.xu] and [self.u1, self.v1] is 1714 constructed such that it is both not close to x and to the boundary 1715 of [self.xl,self.xu] and its function value is estimated 1716 from a local quadratic model around x 1717 i: the index of the point around which the fit is to be computed 1718 1719 The local quadratic model around x is given by 1720 1721 q(y) = self.f[i] + self.gradient*(y-x)' + \ 1722 self.sigma*( (y-x)*diag(D)*(y-x)' + self.df[i] ) 1723 1724 For a row vector y, where D = df0 / dx**2 1725 1726 Output: 1727 ------- 1728 y Point in the intersection of [self.xl,self.xu] and 1729 [self.u1, self.v1] 1730 fy Corresponding estimated function value 1731 """ 1732 # Class 4. 1733 y = vector(self.n) 1734 for j in xrange(self.n): 1735 if self.x[i,j]-self.xl[i,j] > self.xu[i,j]-self.x[i,j]: 1736 y[j] = 0.5*(self.xl[i,j] + self.x[i,j] ) 1737 else: 1738 y[j] = 0.5*(self.x[ i,j] + self.xu[i,j] ) 1739 1740 y = self.round(numpy.minimum( numpy.maximum(y,self.u1),self.v1), 1741 numpy.maximum(self.xl[i,:], self.u1), 1742 numpy.minimum(self.xu[i,:], self.v1) 1743 ) 1744 fy = self._estimate(self.f[i], 1745 self.df[i], 1746 y, 1747 self.x[i], 1748 self.gradient[i,:], 1749 self.sigma[i] 1750 ) 1751 return (y,fy)

1752 1753 1754 #------------------------------------ 1755 # Step 9 1756 #------------------------------------

1757 - def genPointsFromClass234(self):

1758 """ Generate the new points for Class 2,3,4 1759 1760 If X contains any local points, first at most 1761 1762 loc := self.globloc - self.glob 1763 1764 points generated from the local points y are chosen in the order of 1765 ascending model function values fy to yield points of Class 2. 1766 1767 If the desired number of loc points has not been reached yet, 1768 afterwards points y pertaining to nonlocal x which belong to X are 1769 taken (again in the order of ascending fy). 1770 For each potential point of Class 2 or 3 generated in this way, 1771 a subbox [ \underline{x}^j , \bar{x}^j ] of the box tree with 1772 minimal smallness is determined with 1773 y belong to [ \underline{x}^j , \bar{x}^j ] 1774 if 1775 min_i( (\bar{x}^j_i - \underline{x}^j_i )/(v_i-u_i) ) > 1776 0.05 * max_i( (\bar{x}^j_i - \underline{x}^j_i )/(v_i-u_i) ), 1777 1778 the point y is not put into the list of recommended evaluation points 1779 but instead we set self.J4 = self.J4 and {j}, i.e., the box is marked 1780 for the generation of a recommended evaluation point of Class 4. 1781 1782 Note: we denote the X the set of points for which the objective 1783 function has already been evaluated at some stage of SNOBFIT. 1784 """ 1785 loc = self.globloc - self.glob 1786 if loc: 1787 [local,nonlocal] = self._calcLocalPoints( self.ind ) 1788 k = numpy.argsort( self.fy[local] ) 1789 self._genPoints(k, local, classID=2) 1790 1791 if self.request.shape[0] < self.nreq - self.glob: 1792 k = numpy.argsort( self.fy[nonlocal] ) 1793 self._genPoints(k, nonlocal, classID=3) 1794 1795 #--------------------------------------------------- 1796 # There are some suggested evaluation points which don't put into list 1797 # Here we put it into list with Class 4 1798 sreq = self.request.shape[0] 1799 for i in self.J4: 1800 self.ind = removeByInd( self.ind, find(self.ind==i) ) 1801 1802 [y, fy] = self._constructPoints(i) 1803 1804 yy = self.request[:,0:self.n] - \ 1805 numpy.array([y]).repeat(sreq, axis=0) 1806 if numpy.max( abs( y-self.x[i,:])-self.dx) >= -eps and \ 1807 ( not sreq or \ 1808 min( numpy.max(abs(yy) - \ 1809 numpy.array([self.dx]).repeat(sreq, axis=0), 1810 axis=1) ) >= -eps): 1811 1812 self._addRequest( [ self.newRequest(y, fy, 4) ] ) 1813 1814 sreq = self.request.shape[0] 1815 if sreq == self.nreq: 1816 break

1817 1818 #------------------------------------ 1819 # Step 10 1820 #------------------------------------

1821 - def genClass4_FromSmall(self):

1822 """ 1823 Let Smin and Smax denote the minimal and maximal smallness, 1824 respectively, among the boxes in current box tree, and let 1825 1826 M := floor(1/3(Smax-Smin)). 1827 1828 For each smallness S = Smin + m, m = 0, . . . ,M, the boxes are sorted 1829 according to ascending function values f(x) (i.e., self.f( self.x) ). 1830 First a point of Class 4 is generated from the box with S = Smin with 1831 smallest f(x). If self.J4 != Null, then the points of Class 4 1832 belonging to the subboxes with indices in self.J4 are generated. 1833 Subsequently, a point of Class 4 is generated in the box with 1834 smallest f(x) at each nonempty smallness level Smin+m, m = 1, ... ,M, 1835 and then the smallness levels from Smin to Smin+M are gone through 1836 again etc. until we either have nreq recommended evaluation points 1837 or there are no eligible boxes for generating points of Class 4 1838 any more. We assign to the points of Class 4 the model function values 1839 obtained from the local models pertaining to the points in their boxes. 1840 """ 1841 first = 1 1842 sreq = self.request.shape[0] 1843 (Smax,Smin) = self._maxmin( self.small[self.ind] ) 1844 M = int( numpy.floor( (Smax-Smin)/3.0 ) ) 1845 while sreq < self.nreq and len(self.ind) != 0: 1846 for m in xrange(M+1): 1847 1848 if first == 1: 1849 first = 0 1850 continue 1851 1852 m = 0 1853 k = find( self.small[self.ind] == Smin+m ) 1854 while len(k)==0: 1855 m += 1 1856 k = find( self.small[self.ind] == Smin+m ) 1857 1858 if len(k) !=0 : 1859 k = self.ind[k] 1860 k = k[ numpy.argsort(self.f[k]) ] 1861 i = k[0] 1862 self.ind = removeByInd( self.ind, find(self.ind==i) ) 1863 1864 [y, fy] = self._constructPoints(i) 1865 1866 yy = self.request[:,0:self.n] - \ 1867 numpy.array( [y] ).repeat(sreq, axis=0) 1868 dxy = numpy.array([numpy.maximum(self.dy, self.dx)] 1869 ).repeat(sreq, axis=0) 1870 if numpy.max( abs(y-self.x[i,:])-self.dx) >= -eps and \ 1871 ( not sreq or \ 1872 min( numpy.max(abs(yy)-dxy, axis=1) ) >= -eps): 1873 1874 self._addRequest( [self.newRequest(y, fy, 4)] ) 1875 1876 sreq = self.request.shape[0] 1877 if sreq == self.nreq: 1878 break 1879 1880 m = 0

1881 1882 1883 #------------------------------------ 1884 # Step 11 1885 #------------------------------------

1886 - def genPointsFromClass5(self):

1887 """ 1888 If the number of the recommended evaluation points is still 1889 less than self.nreq, the set of evaluation points is filled up 1890 with points of Class 5 as described in Section 4. 1891 If local models are already available, the model function values 1892 for the points of Class 5 are determined as the ones for the points 1893 of Class 4 (see Step 10); otherwise, they are set to NaN. 1894 """ 1895 x1 = self._genClass5Points( 1896 numpy.append( self.x, self.request[:,0:self.n], axis=0 ), 1897 self.nreq - self.request.shape[0] 1898 ) 1899 for j in xrange( x1.shape[0] ): 1900 1901 i=self._sumFind(within(self.xl, 1902 numpy.array([x1[j:]]).repeat(self.nx,axis=0), 1903 self.xu), 1904 self.n 1905 ) 1906 if len(i) > 1: 1907 i = i[ numpy.argmin(self.small[i]) ] 1908 1909 f = self._estimate(self.f[i], 1910 self.df[i], 1911 x1[j,:], 1912 self.x[i,:], 1913 self.gradient[i,:], 1914 self.sigma[i] 1915 ) 1916 1917 self._addRequest([ self.newRequest(x1[j,:],f,5) ])

1918 1919 1920

1921 - def _prepare(self):

1922 """ Some preparation work 1923 1) Init request 1924 2) Make sure that the first a point of Class 4 is generated 1925 from the box with S = Smin with smallest f(x). 1926 etc. 1927 """ 1928 self.request = numpy.zeros( (0, self.n+2) ) 1929 1930 # Make sure that lower and upper bounds of the search boxes 1931 # to belong to the user's specified bounds at begining. 1932 xx = numpy.logical_and( self.xl <= \ 1933 numpy.array([self.v1]).repeat(self.nx,axis=0), 1934 self.xu >= \ 1935 numpy.array([self.u1]).repeat(self.nx,axis=0) ) 1936 self.ind = find( numpy.sum( xx, axis=1 ) == self.n ) 1937 1938 # Find the first index of the recommended eval pts of Class 4 1939 k = find( self.small[self.ind] == min(self.small[self.ind]) ) 1940 k = k[ numpy.argsort(self.f[ self.ind[k] ]) ] 1941 self.J4 = numpy.array( [k[0]] ) 1942 1943 # The curent best solution 1944 [self.fbest, self.jbest] = min_( self.f ) 1945 ind = numpy.argsort(self.f) 1946 self.xbest = self.x[self.jbest,:]

1947 1948

1949 - def _getExtrmeFunc(self):

1950 notnan = find( numpy.isfinite( self.f ) ) 1951 return self._notnanMaxMin(notnan)

1952 1953

1954 - def _isNotEnoughReq(self):

1955 """ 1956 If there is enough request return True. 1957 Otherwise return False 1958 """ 1959 if self.request.shape[0] < self.nreq: 1960 return True 1961 else: 1962 return False

1963 1964

1965 - def _q(self, f, f0, Delta):

1966 return ( f - f0 ) / ( Delta + abs(f - f0) )

1967

1968 - def softmerit(self, f, f0, Delta, x, Fx=None):

1969 """ 1970 Merit function of the soft optimality theorem 1971 1972 Input: 1973 ------ 1974 f objective function value 1975 f0 scalar parameter in the merit function 1976 Delta scalar, positive parameter in the merit function 1977 x the variable. 1978 """ 1979 if Fx is not None: 1980 _Fx = Fx 1981 else: 1982 _Fx = self.constraint.F(x) 1983 1984 if not numpy.isfinite(f) or \ 1985 not numpy.any( numpy.isfinite( _Fx ) ): 1986 # If the objective function or one of the constraint functions is 1987 # infinite or NaN, set the merit function value to 3 1988 return 3 1989 1990 return self._q(f, f0, Delta) + self.constraint.r(x)

1991 1992

1993 - def softmeritOLD(self, f,F,F1,F2,f0,Delta,sigma):

1994 """ 1995 Merit function of the soft optimality theorem 1996 1997 Input: 1998 ------ 1999 f objective function value 2000 F vector containing the values of the constraint functions 2001 (m-vector) 2002 F1 m-vector of lower bounds of the constraints 2003 F2 m-vector of upper bounds of the constraints 2004 f0 scalar parameter in the merit function 2005 Delta scalar, positive parameter in the merit function 2006 sigma positive m-vector, where sigma(i) is the permitted violation 2007 of constraint i 2008 """ 2009 if not numpy.isfinite(f) or \ 2010 not numpy.any( numpy.isfinite(F) ): 2011 # If the objective function or one of the constraint functions is 2012 # infinite or NaN, set the merit function value to 3 2013 fm = 3 2014 return fm 2015 2016 m = len(F) 2017 delta = 0 2018 for i in xrange(m): 2019 if F[i] < F1[i]: 2020 delta += (F1[i]- F[i])**2/sigma[i]**2 2021 2022 elif F[i] > F2[i]: 2023 delta += (F[i] - F2[i])**2/sigma[i]**2 2024 2025 fm = self._q(f, f0, Delta) + 2*delta/(1+delta) 2026 return fm

2027 2028 2029 # -----------------------------------------------------------------

2030 - def _calcCovarAtSolution(self, x, f, df=None):

2031 """ calc the Correlation Matrix or Uncertainty at the solution 2032 2033 OutPut: 2034 ------- 2035 covar the Correlation Matrix. 2036 """ 2037 if x is not None: self.xnew = x.copy() 2038 if f is not None: self.fnew = f.copy() 2039 if df is not None: 2040 self.dfnew = df.copy() 2041 else: 2042 self.dfnew = numpy.ones(len(f))*max(3*self.fac,numpy.sqrt(eps)) 2043 2044 self.updateWorkspace() 2045 2046 # Current best solution 2047 [fSolution, j] = min_( self.f ) 2048 xSolution = self.x[j,:] 2049 2050 K = min( self.nx-1, self.n*(self.n+3) ) 2051 distance=numpy.sum((self.x - \ 2052 numpy.array([xSolution]).repeat(self.nx,axis=0))**2, 2053 axis=1 2054 ) 2055 ind = numpy.argsort(distance)[1:K+1] 2056 2057 # S^k = x^k - xSolution 2058 S = self.x[ind,:] - numpy.array([xSolution]).repeat(K,axis=0) 2059 2060 R = triu( qr(S) ) 2061 L = inv(R).T 2062 2063 # Scale matrix 2064 sc = numpy.sum( numpy.dot(S,L.T)**2, axis=1) ** (3.0/2.0 ) 2065 b = (self.f[ind]-fSolution)/sc 2066 2067 xx = self.x[ind,:] - numpy.array([xSolution]).repeat(K,axis=0) 2068 xx2 = 0.5*xx*xx 2069 A = numpy.bmat( 'xx xx2') 2070 2071 for i in xrange(self.n-1): 2072 xx =( self.x[ind,i] - xSolution[i]) 2073 B = numpy.array( toCol(xx) ).repeat( self.n-i-1, axis=1) 2074 xx = B*(self.x[ind,i+1:self.n] - \ 2075 numpy.array([ xSolution[i+1:self.n] ]).repeat(K,axis=0) ) 2076 A = numpy.bmat('A xx') 2077 2078 xx = numpy.array( toCol(sc) ).repeat(A.shape[1], axis=1) 2079 A = A/xx 2080 2081 # In q, it includes the components of g and the upper triangle of G 2082 q = mldiv(A, toCol(b)) 2083 2084 # Jacobian matrix 2085 J = numpy.diag( q[self.n:2*self.n] )/2 2086 k = 2*self.n 2087 for i in xrange( self.n-1 ): 2088 for j in xrange(self.n-i-1): 2089 J[i, j+i+1] = q[k]/2 2090 J[j+i+1, i] = q[k]/2 2091 k += 1 2092 2093 # The covariance matrix is computed by, 2094 # covar = (J^T J)^{-1} 2095 covar = inv(J.T * J) 2096 2097 # If the minimisation uses the weighted least-squares function. 2098 # the covariance matrix should be multiplied by the variance of 2099 # the residuals about the best-solution least-squares. 2100 # For example, the fit of the reflectometry problem. 2101 if self.isLeastSqrt: 2102 # Calc the squared residuals(i.e., the goodness of fit) 2103 # at the solution. 2104 sumsq = abs(fSolution)**2/self.n 2105 covar *= sumsq 2106 2107 2108 return covar

2109 2110

2111 - def _covarianceMatrix(self):

2112 """ 2113 calculate the Covariance Matrix at the solution 2114 """ 2115 # Compute function values at the suggested points 2116 nrequest = self.request.shape[0] 2117 x = numpy.zeros( (nrequest,self.n) ) 2118 f = vector( nrequest ) 2119 for j in xrange( nrequest ): 2120 x[j,:] = self.request[j,0:self.n] 2121 f[j] = self.func(x[j,:]) 2122 2123 return self. _calcCovarAtSolution(x, f)

2124 2125 # The alias 2126 covar = _covarianceMatrix 2127 2128

2129 - def uncertainty(self):

2130 """ 2131 calculate the uncertainty of the fitting parameter at the solution. 2132 """ 2133 covar = self._covarianceMatrix() 2134 2135 # calculate the uncertaity of fit parameters. 2136 uct = numpy.sqrt( numpy.diag(covar) ) 2137 2138 return uct

2139 2140

2141 - def fit(self, x=None, f=None, df=None, initialCall=False):

2142 """ 2143 One call of snobfit 2144 2145 Output: 2146 ------- 2147 request nreq x (n+3)-matrix 2148 request[j,0:n] is the jth newly generated point, 2149 request[j,n+1] is its estimated function value and 2150 request[j,n+2] is its uncertainty of estimated function value 2151 request[j,n+3] indicates for which reason the point 2152 = 1 best prediction 2153 = 2 putative local minimizer 2154 = 3 alternative good point 2155 = 4 explore empty region 2156 = 5 to fill up the required number of function values 2157 if too little points of other classes are found 2158 xbest current best point 2159 fbest current best function value (i.e. function value at xbest) 2160 """ 2161 if initialCall: 2162 # When a new run is started 2163 self.update_uv() # Step 1 2164 self.update_input() # Step 2 2165 self.branch() # step 3 2166 2167 (fmax, fmin) = self._getExtrmeFunc() 2168 if self.nx >= self.snn + 1 and fmin < fmax: 2169 self.calc_snn() # Step 4 2170 self.calc_fdf_nan() # step 5 2171 self.localFit() # Step 6 2172 else: 2173 # This is equal go to Step 11 2174 return self.returnNanRequest() 2175 2176 else: # a continue call 2177 # For the continue call, we needn't go through Step 1-6, 2178 # We just update the work space for last call 2179 if x is not None: self.xnew = x.copy() 2180 if f is not None: self.fnew = f.copy() 2181 if df is not None: 2182 self.dfnew = df.copy() 2183 else: 2184 self.dfnew = numpy.ones(len(f))*max(3*self.fac,numpy.sqrt(eps)) 2185 2186 self.updateWorkspace() 2187 2188 if isEmpty(self.near): # This is equal go to Step 11 2189 self.returnNanRequest() 2190 2191 self._prepare() 2192 2193 self.locQuadFit() # Step 7 2194 2195 if self._isNotEnoughReq(): 2196 self.assignN4EachClass() # Step 8 2197 self.genPointsFromClass234() # step 9 2198 2199 self.genClass4_FromSmall() # step 10 2200 2201 if self._isNotEnoughReq(): 2202 self.genPointsFromClass5() # Step 11 2203 2204 if self._isNotEnoughReq(): 2205 self._warning()

2206 2207 2208

2209 - def solve(self):

2210 """ Main Loop 2211 2212 The optimization is stopped at one of the follow conditions: 2213 1: if approximately ncall fucntion values have been exceeded. 2214 2215 2: if at least minfval function values were obtained and 2216 the best function value wasn't improved in the last nstop 2217 calls to SNOBFIT 2218 2219 3: stop when within tolerance of the known global minimum. E.g., 2220 for least squares problems, the global minimum is 0 for a 2221 perfect match. Various test functions also have known global 2222 minima. 2223 2224 Output: 2225 ------- 2226 xopt minimizer of function. 2227 fopt value of function at minimum: fopt = func(xopt). 2228 ncall number of iterations. 2229 """ 2230 # Function call counter 2231 ncall = self.nreq 2232 xopt = inf 2233 improvement = 0 2234 2235 # iteration counter 2236 i = 0 2237 2238 # Define x and f 2239 x = numpy.zeros( (self.nreq, self.n) ) 2240 f = numpy.zeros( self.nreq ) 2241 2242 # Repeat until the limit on function calls is reached 2243 while ncall < self.maxiter: 2244 if ncall == self.nreq: 2245 # initial call 2246 self.fit(initialCall=True) 2247 2248 else: 2249 # continuation call 2250 self.fit(x=x, f=f) 2251 2252 [xopt, fopt] = (self.xbest, self.fbest) 2253 # Compute (perturbed) function values at the suggested points 2254 for j in xrange( self.nreq ): 2255 x[j,:] = self.request[j,0:self.n] 2256 f[j] = self.func(x[j,:]) 2257 2258 # update function call counter 2259 ncall += self.nreq 2260 2261 # best of the new function values 2262 [fbestn, jbest] = min_(f) 2263 2264 # If a better function value has been found, update fbest 2265 if fbestn < fopt: 2266 fopt = fbestn 2267 xopt = x[jbest,:] 2268 improvement = 0 2269 else: 2270 improvement += 1 2271 2272 # Stop if no improvement expected 2273 if improvement >= self.nstop: 2274 break 2275 2276 if self.callback is not None : 2277 self.callback(i, x[jbest], fbestn, improvement==0) 2278 2279 # Stop at user specified stopping criterion 2280 # best function value is found with a relative error < 1.e-6 2281 if self.fglob is not None: 2282 if self.fglob: 2283 if abs((fopt-self.fglob)/self.fglob) < self.rtol: 2284 if self.disp: print "The best solution:", xopt,fopt 2285 break 2286 else: # safeguard if functions with fglob=0 are added by user 2287 if abs(fopt) < self.ftol: 2288 if self.disp: print "The best results", xopt, fopt 2289 break 2290 2291 if self.retall: 2292 print xopt, fopt, ncall 2293 2294 i+=1 2295 2296 if self.retuct: 2297 return xopt, self.uncertainty(), fopt, ncall 2298 else: 2299 return xopt, fopt, ncall

2300 2301

2302 - def softSolve(self):

2303 """ Main Loop 2304 2305 The optimization is stopped at one of the follow conditions: 2306 1: if approximately ncall fucntion values have been exceeded. 2307 2308 2: if at least minfval function values were obtained and 2309 the best function value wasn't improved in the last nstop 2310 calls to SNOBFIT 2311 2312 3: stop when within tolerance of the known global minimum. E.g., 2313 for least squares problems, the global minimum is 0 for a 2314 perfect match. Various test functions also have known global 2315 minima. 2316 2317 Output: 2318 ------- 2319 xopt minimizer of function. 2320 fopt value of function at minimum: fopt = func(xopt). 2321 ncall number of iterations. 2322 """ 2323 # Function call counter 2324 ncall = self.nreq 2325 xopt = inf 2326 improvement = 0 2327 change = 0 2328 2329 # iteration counter 2330 i = 0 2331 2332 # Calc the objective function for the starting points. 2333 _x = self.x.copy() 2334 _f = vector( self.nreq ) 2335 for j in xrange( self.nreq ): 2336 _f[j] = self.func( _x[j] ) 2337 2338 2339 # Calc the constraints for the starting points. 2340 # get scalar parameter in the merit function 2341 _F = numpy.zeros( ( self.nreq, len(self.constraint.C) ) ) 2342 f0 = self._f0 2343 Delta = self._Delta 2344 for j in xrange( len(self.f) ): 2345 _F[j] = self.constraint.F(self.x[j,:]) 2346 2347 # Some declaration of vector 2348 F = numpy.zeros( (self.nreq, len(self.constraint.C)) ) 2349 f = numpy.zeros( self.nreq ) 2350 2351 # Two common constants 2352 sigmaLo = self.constraint.sigmaLo() 2353 sigmaHi = self.constraint.sigmaHi() 2354 2355 # Repeat until the limit on function calls is reached 2356 while ncall < self.maxiter: 2357 if ncall == self.nreq: 2358 # initial call 2359 self.fit(initialCall=True) 2360 2361 else: 2362 # continuation call 2363 self.fit(x=x, f=fm) 2364 2365 [xopt, fopt] = (self.xbest, self.fbest) 2366 2367 # Compute (perturbed) function values at the suggested points 2368 x = numpy.zeros( (self.nreq, self.n) ) 2369 fm = numpy.zeros( self.nreq ) 2370 2371 for j in xrange( self.nreq ): 2372 x[j,:] = self.request[j,0:self.n] 2373 f[j] = self.func(x[j,:]) 2374 FF = self.constraint.F(x[j,:]) 2375 F[j] = FF 2376 fm[j] = self.softmerit(f[j] , f0, Delta, x[j,:], Fx=FF) 2377 2378 if f[j] <= self.fglob and \ 2379 min(FF - self.F1 + sigmaLo) >= 0 and \ 2380 min(self.F2 + sigmaHi -FF) >= 0: 2381 ncall += j 2382 xsoft = x[j,:] 2383 [fbestn,jbest] = min_(fm) 2384 if fbestn < fopt: 2385 fopt = fbestn 2386 xopt = x[jbest,:] 2387 2388 # show number of function values, best point and function 2389 # value and the point xsoft fulfilling the conditions of 2390 # the soft optimality theorem (hopefully close to xglob) 2391 if self.disp: 2392 print ncall,xopt,fopt,xsoft 2393 2394 # Append new suggested points 2395 _x = numpy.append(_x, x, axis=0) 2396 _f = numpy.append(_f, f, axis=0) 2397 _F = numpy.append(_F, F, axis=0) 2398 2399 # update function call counter 2400 ncall += self.nreq 2401 2402 # best of the new function values 2403 [fbestn, jbest] = min_(fm) 2404 2405 # If a better function value has been found, update fbest 2406 if fbestn < fopt: 2407 fopt = fbestn 2408 xopt = x[jbest,:] 2409 improvement = 0 2410 else: 2411 improvement += 1 2412 2413 # Stop if no improvement expected 2414 if improvement >= self.nstop: 2415 break 2416 2417 if self.callback is not None : 2418 self.callback(i, x[jbest], fbestn, improvement==0) 2419 2420 if fopt < 0 and change == 0: 2421 n = _x.shape[0] 2422 c1 = numpy.min( _F - numpy.array([self.F1]).repeat(n,axis=0), 2423 axis=1 ) > -eps 2424 c2 = numpy.min( numpy.array([self.F2]).repeat(n,axis=0)-_F, 2425 axis=1) > -eps 2426 ind = find( numpy.logical_and(c1, c2) ) 2427 if len(ind) != 0 : 2428 change = 1 2429 f0 = min(_f[ind]) 2430 Delta = numpy.median( abs(f-f0) ) 2431 fm = vector( n ) 2432 for j in xrange(n): 2433 fm[j]=self.softmerit(_f[j], f0, Delta, _x[j] ) 2434 2435 x = _x.copy() 2436 2437 _fopt = self.func(xopt) 2438 2439 # Stop at user specified stopping criterion 2440 # best function value is found with a relative error < 1.e-6 2441 if self.fglob is not None: 2442 if self.fglob: 2443 if abs((_fopt-self.fglob)/self.fglob) < self.rtol: 2444 if self.disp: print "The best solution:", xopt,_fopt 2445 break 2446 else: # safeguard if functions with fglob=0 are added by user 2447 if abs(_fopt) < self.ftol: 2448 if self.disp: print "The best results", xopt, _fopt 2449 break 2450 2451 if self.retall: 2452 print xopt, _fopt, ncall 2453 2454 i+=1 2455 2456 if self.retuct: 2457 return xopt, self.uncertainty(), _fopt, ncall 2458 else: 2459 return xopt, _fopt, ncall

Source Code for Module snobfit.snobfit