Changeset 171
- Timestamp:
- 08/07/09 11:27:13 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/alta/mystic-0.1a2/mystic/snobfit_solver.py
r168 r171 27 27 - StepMonitor = Sow() 28 28 - enable_signal_handler() 29 - termination = SnobfitTermination(rtol, ftol, nstop)29 - termination = SnobfitTermination(rtol, ftol, generations) 30 30 31 31 Usage … … 72 72 # Import Mystic tools 73 73 from mystic.tools import Null, wrap_function 74 from mystic.tools import wrap_bounds , random_seed74 from mystic.tools import wrap_bounds 75 75 76 76 from abstract_solver import AbstractSolver … … 123 123 p -- probability of generating a evaluation point of Class 4. 124 124 fac -- Factor for multiplicative perturbation of the data. 125 disp -- non-zero if fval and warnflag outputs are desired.126 callback -- an optional user-supplied function to call after each127 iteration. It's called as callback(xbest), where xbest128 is the current best point.129 125 fglob -- the user specified global function value. 130 126 … … 134 130 self.p = 0.5 135 131 self.fac = 0.0 136 disp = 0137 callback = None138 132 self.fglob = 0 139 133 … … 141 135 if kwds.has_key('p'): self.p = kwds['p'] 142 136 if kwds.has_key('fac'): self.fac = kwds['fac'] 143 if kwds.has_key('disp'): disp = kwds['disp']144 if kwds.has_key('callback'): callback = kwds['callback']145 137 if kwds.has_key('fglob'): self.fglob = kwds['fglob'] 146 138 … … 160 152 #------------------------------------------------------------- 161 153 154 self.func = func 155 162 156 x0 = self.population[0] 163 157 self.x = self.population … … 182 176 183 177 # Function call counter 184 ncall = self.nPop 178 #ncall = self.nPop # <-- original snobfit's function call counter 179 initfcalls = fcalls[0] 185 180 186 181 # iteration counter … … 195 190 196 191 # Repeat until the limit on function calls is reached 197 while ncall < self._maxiter: 192 #while ncall < self._maxiter: 193 while fcalls[0] < self._maxiter: 198 194 StepMonitor(self.bestSolution, self.bestEnergy) 199 200 if ncall == self.nPop: 195 196 #if ncall == self.nPop: 197 if fcalls[0] == initfcalls: 201 198 # initial call 202 199 self.fit(func, initialCall=True) … … 213 210 214 211 # update function call counter 215 ncall += self.nPop212 #ncall += self.nPop 216 213 217 214 # best of the new function values … … 223 220 xopt = x[jbest,:] 224 221 225 if callback is not None:226 callback(x[jbest])222 if self.callback is not None: 223 self.callback(x[jbest]) 227 224 228 225 # Monitor progress … … 245 242 246 243 if fcalls[0] >= self._maxfun: 247 if disp:244 if self.disp: 248 245 print "Warning: Maximum number of function evaluations has "\ 249 246 "been exceeded." 250 247 elif i >= self._maxiter: 251 if disp:248 if self.disp: 252 249 print "Warning: Maximum number of iterations has been exceeded" 253 250 else: 254 if disp:251 if self.disp: 255 252 print "Optimization terminated successfully." 256 253 print " Current function value: %f" % fval … … 259 256 return 260 257 261 # Overwrite abstract_solver's SetInitialPoints() 262 def SetInitialPoints(self, x0): 263 """Set initial points with guess x0.""" 264 self.population = self._setInitRecommendedEvalPoints(x0) 265 258 # Overwrite abstract_solver's SetRandomInitialPoints() 259 def SetRandomInitialPoints(self, min=None, max=None): 260 """Set random initial points.""" 261 self.population = self._setInitRecommendedEvalPoints() 262 263 # def SetInitialPoints(self, x0): 264 # self.SetRandomInitialPoints() 265 # self.population[0] = x0 266 267 266 268 ################################################################################# 267 269 268 270 # Helper methods for SnobfitSolver 269 271 270 def _setInitRecommendedEvalPoints(self , x0):272 def _setInitRecommendedEvalPoints(self): 271 273 """ 272 274 From the initial user-supplied guess, we set up the first set of … … 277 279 Make sure the initial guess x0 is in recommended evaluation points. 278 280 """ 279 x = numpy.dot( rand(self.nPop -1, self.nDim),281 x = numpy.dot( rand(self.nPop, self.nDim), 280 282 numpy.diag(self._strictMax-self._strictMin) 281 283 ) + \ 282 numpy.array( [self._strictMin] ).repeat(self.nPop -1,axis=0)283 284 x = numpy.append( [x0], x, axis=0 )284 numpy.array( [self._strictMin] ).repeat(self.nPop,axis=0) 285 286 #x = numpy.append( [x0], x, axis=0 ) 285 287 return x 286 288 … … 2163 2165 self._warning() 2164 2166 2167 ########################################################################### 2168 # Helper methods for SnobfitSolver to calculate the covariance 2169 # matrix, uncertainty: 2170 2171 def _calcCovarAtSolution(self, x, f, df=None): 2172 """ calc the Correlation Matrix or Uncertainty at the solution 2173 2174 OutPut: 2175 ------- 2176 covar the Correlation Matrix. 2177 """ 2178 if x is not None: self.xnew = x.copy() 2179 if f is not None: self.fnew = f.copy() 2180 if df is not None: 2181 self.dfnew = df.copy() 2182 else: 2183 self.dfnew = numpy.ones(len(f))*max(3*self.fac,numpy.sqrt(eps)) 2184 2185 self.updateWorkspace(self.func) 2186 2187 # Current best solution 2188 [fSolution, j] = min_( self.f ) 2189 xSolution = self.x[j,:] 2190 2191 K = min( self.nx-1, self.nDim*(self.nDim+3) ) 2192 distance=numpy.sum((self.x - \ 2193 numpy.array([xSolution]).repeat(self.nx,axis=0))**2, 2194 axis=1 2195 ) 2196 ind = numpy.argsort(distance)[1:K+1] 2197 2198 # S^k = x^k - xSolution 2199 S = self.x[ind,:] - numpy.array([xSolution]).repeat(K,axis=0) 2200 2201 #---------- Fixed typo ------------------- 2202 #R = numpy.linalg.qr(a, mode='r') 2203 R = numpy.linalg.qr(S, mode='r') 2204 L = inv(R).T 2205 2206 # Scale matrix 2207 sc = numpy.sum( numpy.dot(S,L.T)**2, axis=1) ** (3.0/2.0 ) 2208 b = (self.f[ind]-fSolution)/sc 2209 2210 xx = self.x[ind,:] - numpy.array([xSolution]).repeat(K,axis=0) 2211 xx2 = 0.5*xx*xx 2212 A = numpy.bmat( 'xx xx2') 2213 2214 for i in xrange(self.nDim-1): 2215 xx =( self.x[ind,i] - xSolution[i]) 2216 B = numpy.array( toCol(xx) ).repeat( self.nDim-i-1, axis=1) 2217 xx = B*(self.x[ind,i+1:self.nDim] - \ 2218 numpy.array([ xSolution[i+1:self.nDim] ]).repeat(K,axis=0) ) 2219 A = numpy.bmat('A xx') 2220 2221 xx = numpy.array( toCol(sc) ).repeat(A.shape[1], axis=1) 2222 A = A/xx 2223 2224 # In q, it includes the components of g and the upper triangle of G 2225 q = mldiv(A, toCol(b)) 2226 2227 # Jacobian matrix 2228 J = numpy.diag( q[self.nDim:2*self.nDim] )/2 2229 k = 2*self.nDim 2230 for i in xrange( self.nDim-1 ): 2231 for j in xrange(self.nDim-i-1): 2232 J[i, j+i+1] = q[k]/2 2233 J[j+i+1, i] = q[k]/2 2234 k += 1 2235 2236 # The covariance matrix is computed by, 2237 # covar = (J^T J)^{-1} 2238 covar = inv(J.T * J) 2239 2240 # If the minimisation uses the weighted least-squares function. 2241 # the covariance matrix should be multiplied by the variance of 2242 # the residuals about the best-solution least-squares. 2243 # For example, the fit of the reflectometry problem. 2244 if self.isLeastSqrt: 2245 # Calc the squared residuals(i.e., the goodness of fit) 2246 # at the solution. 2247 sumsq = abs(fSolution)**2/self.nDim 2248 covar *= sumsq 2249 2250 2251 return covar 2252 2253 2254 def _covarianceMatrix(self): 2255 """ 2256 Calculate the Covariance Matrix at the solution 2257 """ 2258 # Compute function values at the suggested points 2259 nrequest = self.request.shape[0] 2260 x = numpy.zeros( (nrequest,self.nDim) ) 2261 f = vector( nrequest ) 2262 for j in xrange( nrequest ): 2263 x[j,:] = self.request[j,0:self.nDim] 2264 f[j] = self.func(x[j,:]) 2265 2266 return self. _calcCovarAtSolution(x, f) 2267 2268 # The alias 2269 covar = _covarianceMatrix 2270 2271 2272 def uncertainty(self, isLeastSqrt=False): 2273 """ 2274 Calculate the uncertainty of the fitting parameter at the solution. 2275 2276 Optional input: isLeastSqrt -- boolean -- True if the minimisation 2277 uses the least squares function. 2278 """ 2279 self.isLeastSqrt=isLeastSqrt 2280 covar = self._covarianceMatrix() 2281 print covar 2282 print numpy.diag(covar) 2283 2284 # calculate the uncertainty of fit parameters. 2285 uct = numpy.sqrt( numpy.diag(covar) ) 2286 2287 return uct 2288 2289 2290 2165 2291 ############################################################################### 2166 # If using an instance of SoftConstraints, use SnobfitSoftSolver instead of SnobfitSolver 2292 # If using an instance of SoftConstraints (constraint is not None), 2293 # use SnobfitSoftSolver instead of SnobfitSolver 2167 2294 2168 2295 class SnobfitSoftSolver(SnobfitSolver): 2169 """Inherits from SnobfitSolver. Use if there are soft constraints.""" 2296 """Inherits from SnobfitSolver. Use this solver only if 2297 there are soft constraints.""" 2170 2298 2171 2299 def Solve(self, func, termination, sigint_callback=None, … … 2191 2319 p -- probability of generating a evaluation point of Class 4. 2192 2320 fac -- Factor for multiplicative perturbation of the data. 2193 disp -- non-zero if fval and warnflag outputs are desired.2194 callback -- an optional user-supplied function to call after each2195 iteration. It's called as callback(xbest), where xbest2196 is the current best point.2197 2321 fglob -- the user specified global function value. 2198 2322 … … 2202 2326 self.p = 0.5 2203 2327 self.fac = 0.0 2204 disp = 02205 callback = None2206 2328 self.fglob = 0 2207 2329 … … 2209 2331 if kwds.has_key('p'): self.p = kwds['p'] 2210 2332 if kwds.has_key('fac'): self.fac = kwds['fac'] 2211 if kwds.has_key(' disp'): disp = kwds['disp']2212 if kwds.has_key('callback'): callback = kwds['callback'] 2333 if kwds.has_key('fglob'): self.fglob = kwds['fglob'] 2334 2213 2335 # SoftConstraints requires a constraint passed to Solve().... 2214 2336 if kwds.has_key('constraint'): self.constraint = kwds['constraint'] 2215 if kwds.has_key('fglob'): self.fglob = kwds['fglob']2216 2337 2217 2338 #------------------------------------------------------------- … … 2253 2374 ncall = self.nPop 2254 2375 xopt = inf 2255 improvement = 02256 2376 change = 0 2257 2377 … … 2342 2462 xopt = x[jbest,:] 2343 2463 2344 if callback is not None :2345 callback(i, x[jbest], fbestn, improvement==0)2464 if self.callback is not None : 2465 self.callback(x[jbest]) 2346 2466 2347 2467 if fopt < 0 and change == 0: … … 2399 2519 print " Function evaluations: %d" % fcalls[0] 2400 2520 return 2401 2402 # Overwrite abstractsolver's SetInitialPoints()2403 def SetInitialPoints(self, x0):2404 """Set initial points with guess x0."""2405 self.population = self._setInitRecommendedEvalPoints(x0)2406 2521 2407 2522 ####################################################################### … … 2488 2603 retall = 0, 2489 2604 callback = None, 2490 seed = None,2491 2605 constraint = None, 2492 full_output = 0): 2606 full_output = 0, 2607 retuct = False): 2493 2608 """Minimize a function using Snobfit - Stable Noisy Optimization by Branch and Fit. 2494 2609 … … 2519 2634 iteration. It's called as callback(xbest), where xbest 2520 2635 is the current best point. 2521 seed -- The seed for numpy.random. If the user set this number,it makes2522 sure that it returns the same solution from repeated tests2523 2636 constraint -- an instance of SoftConstraints. If constraint is not None, then use 2524 2637 SnobfitSoftSolver instead of SnobfitSolver. 2525 2638 full_output -- non-zero if fval and warnflag outputs are 2526 2639 desired. 2527 2528 2529 Returns: (xopt, {fopt, iter, funcalls, warnflag}, {allvecs}) 2640 retuct -- True to return the uncertainty of the fitting parameters 2641 2642 2643 Returns: (xopt, {fopt, iter, funcalls, warnflag}, {allvecs}, {uncertainty}) 2530 2644 2531 2645 xopt -- ndarray - minimizer of function … … 2537 2651 2 : 'Maximum number of iterations.' 2538 2652 allvecs -- list - a list of solutions at each iteration 2653 uncertainty -- list - Returned if retuct is True. The uncertainty of 2654 the fitting parameters. 2539 2655 2540 2656 """ … … 2545 2661 from mystic.termination import SnobfitTermination 2546 2662 2547 # Set the random seed2548 if seed:2549 random_seed(seed)2550 2551 # Messy?:2552 # Using npop instead of dn may not be the best idea?2553 2663 if constraint is None: 2554 2664 if npop: … … 2562 2672 solver = SnobfitSoftSolver(len(x0)) 2563 2673 2674 solver.disp = disp 2675 solver.callback = callback 2676 2564 2677 # Problem?: SetStrictRanges must be called before SetInitialPoints 2565 # In fact, SetStrictRanges must always be called. Perhaps2566 # integrate bounds into __init__? 2678 # In fact, _strictMin and _strictMax must always be set. Perhaps 2679 # integrate bounds into __init__? Make a default? 2567 2680 solver.SetStrictRanges(bounds[0],bounds[1]) 2568 2681 solver.SetInitialPoints(x0) … … 2572 2685 EvaluationMonitor=evalmon,StepMonitor=stepmon,\ 2573 2686 ExtraArgs=args,\ 2574 p=p, fac=fac, disp=disp, fglob=fglob,\ 2575 callback=callback, 2576 constraint=constraint) 2687 p=p, fac=fac, fglob=fglob,constraint=constraint) 2577 2688 solution = solver.Solution() 2578 2689 2579 # code below here pushes output to scipy.optimize interface2580 2690 x = solver.bestSolution 2581 2691 fval = solver.bestEnergy … … 2601 2711 retlist = (x, allvecs) 2602 2712 2713 if retuct: 2714 retlist += (solver.uncertainty(),) 2715 2603 2716 return retlist 2604 2717
Note: See TracChangeset
for help on using the changeset viewer.