root/trunk/openopt/scikits/openopt/Kernel/Point.py

Revision 1839, 19.8 kB (checked in by dmitrey.kroshko, 2 years ago)

fix for recursive connectOOVars() + some other changes

Line 
1 # created by Dmitrey
2 from numpy import copy, isnan, array, argmax, abs, zeros, any, isfinite, all, where, asscalar, sign, dot, sqrt, array_equal
3 __docformat__ = "restructuredtext en"
4 empty_arr = array(())
5
6 class Point:
7     """
8     the class is used to prevent calling non-linear constraints more than once
9     f, c, h are funcs for obtaining objFunc, non-lin ineq and eq constraints.
10     df, dc, dh are funcs for obtaining 1st derivatives.
11     """
12     __expectedArgs__ = ['x', 'f', 'mr']
13     def __init__(self, p, x, *args, **kwargs):
14         self.p = p
15         self.x = copy(x)
16         for i, arg in enumerate(args):
17             setattr(self, '_' + self.__expectedArgs__[i], args[i])
18         for name, val in kwargs.iteritems():
19             setattr(self, '_' + name, val)
20         #assert self.x is not None
21
22     def f(self):
23         if not hasattr(self, '_f'): self._f = self.p.f(self.x)
24         return copy(self._f)
25
26     def df(self):
27         if not hasattr(self, '_df'): self._df = self.p.df(self.x)
28         return copy(self._df)
29
30     def c(self, ind=None):
31         if not self.p.userProvided.c: return empty_arr.copy()
32         if ind is None:
33             if not hasattr(self, '_c'): self._c = self.p.c(self.x)
34             return copy(self._c)
35         else:
36             if hasattr(self, '_c'): return copy(self._c[ind])
37             else: return copy(self.p.c(self.x, ind))
38
39
40     def dc(self, ind=None):
41         if not self.p.userProvided.c: return empty_arr.copy().reshape(0, self.p.n)
42         if ind is None:
43             if not hasattr(self, '_dc'): self._dc = self.p.dc(self.x)
44             return copy(self._dc)
45         else:
46             if hasattr(self, '_dc'): return copy(self._dc[ind])
47             else: return copy(self.p.dc(self.x, ind))
48
49
50     def h(self, ind=None):
51         if not self.p.userProvided.h: return empty_arr.copy()
52         if ind is None:
53             if not hasattr(self, '_h'): self._h = self.p.h(self.x)
54             return copy(self._h)
55         else:
56             if hasattr(self, '_h'): return copy(self._h[ind])
57             else: return copy(self.p.h(self.x, ind))
58
59     def dh(self, ind=None):
60         if not self.p.userProvided.h: return empty_arr.copy().reshape(0, self.p.n)
61         if ind is None:
62             if not hasattr(self, '_dh'): self._dh = self.p.dh(self.x)
63             return copy(self._dh)
64         else:
65             if hasattr(self, '_dh'): return copy(self._dh[ind])
66             else: return copy(self.p.dh(self.x, ind))
67
68     def d2f(self):
69         if not hasattr(self, '_d2f'): self._d2f = self.p.d2f(self.x)
70         return copy(self._d2f)
71
72     def lin_ineq(self):
73         if not hasattr(self, '_lin_ineq'): self._lin_ineq = self.p.__get_AX_Less_B_Residuals__(self.x)
74         return copy(self._lin_ineq)
75
76     def lin_eq(self):
77         if not hasattr(self, '_lin_eq'): self._lin_eq = self.p.__get_AeqX_eq_Beq_Residuals__(self.x)
78         return copy(self._lin_eq)
79
80     def __all_lin_ineq(self):
81         if not hasattr(self, '_all_lin_ineq'):
82             lb, ub, lin_ineq = self.lb(), self.ub(), self.lin_ineq()
83             r = 0
84             # TODO: CHECK IT - when 0 (if some nans), when contol
85             threshold = 0
86 #            if all(isfinite(self.f())): threshold = self.p.contol
87 #            else: threshold = 0
88
89             lb, ub = self.lb(), self.ub()
90             lin_ineq = self.lin_ineq()
91             lin_eq = self.lin_eq()
92             ind_lb, ind_ub = where(lb>threshold)[0], where(ub>threshold)[0]
93             ind_lin_ineq = where(lin_ineq>threshold)[0]
94             ind_lin_eq = where(abs(lin_eq)>threshold)[0]
95
96             if ind_lb.size != 0:
97                 r += sum(lb[ind_lb] ** 2)
98             if ind_ub.size != 0:
99                 r += sum(ub[ind_ub] ** 2)
100             if ind_lin_ineq.size != 0:
101                 r += sum(lin_ineq[ind_lin_ineq] ** 2)
102             if ind_lin_eq.size != 0:
103                 r += sum(lin_eq[ind_lin_eq] ** 2)
104             self._all_lin_ineq = sqrt(r)
105         return copy(self._all_lin_ineq)
106
107     def __all_lin_ineq_gradient(self):
108         if not hasattr(self, '_all_lin_ineq_gradient'):
109             p = self.p
110             d = zeros(p.n)
111             contol = p.contol
112
113             # TODO: CHECK IT - when 0 (if some nans), when contol
114 #            if all(isfinite(self.f())): threshold = contol
115 #            else: threshold = 0
116             threshold = 0
117
118             lb, ub = self.lb(), self.ub()
119             lin_ineq = self.lin_ineq()
120             lin_eq = self.lin_eq()
121             ind_lb, ind_ub = where(lb>threshold)[0], where(ub>threshold)[0]
122             ind_lin_ineq = where(lin_ineq>threshold)[0]
123             ind_lin_eq = where(abs(lin_eq)>threshold)[0]
124
125             if ind_lb.size != 0:
126                 d[ind_lb] -= lb[ind_lb]# d/dx((x-lb)^2) for violated constraints
127             if ind_ub.size != 0:
128                 d[ind_ub] += ub[ind_ub]# d/dx((x-ub)^2) for violated constraints
129             if ind_lin_ineq.size != 0:
130                 a = p.A[ind_lin_ineq]
131                 b = p.b[ind_lin_ineq]
132                 d += dot(a.T, dot(a, self.x)  - b) # d/dx((Ax-b)^2)
133             if ind_lin_eq.size != 0:
134                 aeq = p.Aeq[ind_lin_eq]
135                 beq = p.beq[ind_lin_eq]
136                 d += dot(aeq.T, dot(aeq, self.x)  - beq) # 0.5*d/dx((Aeq x - beq)^2)
137             devider = self.__all_lin_ineq()
138             if devider != 0:
139                 self._all_lin_ineq_gradient = d / devider
140             else:
141                 self._all_lin_ineq_gradient = d
142         return copy(self._all_lin_ineq_gradient)
143
144     def lb(self):
145         if not hasattr(self, '_lb'): self._lb = self.p.lb - self.x
146         return copy(self._lb)
147
148     def ub(self):
149         if not hasattr(self, '_ub'): self._ub = self.x - self.p.ub
150         return copy(self._ub)
151
152     def mr(self, retAll = False):
153         # returns max residual
154         return self.__mr(retAll)
155
156     def __mr(self, retAll = False):
157         if not hasattr(self, '_mr'):
158             r, fname, ind = 0, None, None
159             for field in ('c',  'lin_ineq', 'lb', 'ub'):
160                 fv = array(getattr(self, field)()).flatten()
161                 if fv.size > 0:
162                     ind_max = argmax(fv)
163                     val_max = fv[ind_max]
164                     if r < val_max:
165                         r, ind, fname = val_max, ind_max, field
166             for field in ('h', 'lin_eq'):
167                 fv = array(getattr(self, field)()).flatten()
168                 if fv.size > 0:
169                     fv = abs(fv)
170                     ind_max = argmax(fv)
171                     val_max = fv[ind_max]
172                     if r < val_max:
173                         r, ind, fname = val_max, ind_max, field
174             self._mr, self._mrName,  self._mrInd= r, fname, ind
175         if retAll:
176             return asscalar(copy(self._mr)), self._mrName, asscalar(copy(self._mrInd))
177         else: return asscalar(copy(self._mr))
178
179     def mr_alt(self, retAll = False):
180         if not hasattr(self, '_mr_alt'):
181             mr, fname, ind = self.__mr(retAll = True)
182             self._mr_alt, self._mrName_alt,  self._mrInd_alt= mr, fname, ind
183             c, h= self.c(), self.h()
184             all_lin_ineq = self.__all_lin_ineq()
185             r = 0
186             Type = 'all_lin_ineq'
187             if c.size != 0:
188                 ind_max = argmax(c)
189                 val_max = c[ind_max]
190                 if val_max > r:
191                     r = val_max
192                     Type = 'c'
193                     ind = ind_max
194             if h.size != 0:
195                 h = abs(h)
196                 ind_max = argmax(h)
197                 val_max = h[ind_max]
198                 #hm = abs(h).max()
199                 if val_max > r:
200                     r = val_max
201                     Type = 'h'
202                     ind = ind_max
203 #            if lin_eq.size != 0:
204 #                l_eq = abs(lin_eq)
205 #                ind_max = argmax(l_eq)
206 #                val_max = l_eq[ind_max]
207 #                if val_max > r:
208 #                    r = val_max
209 #                    Type = 'lin_eq'
210 #                    ind = ind_max
211
212             if  r <= all_lin_ineq:
213                 self._mr_alt, self._mrName_alt,  self._mrInd_alt = all_lin_ineq, 'all_lin_ineq', 0
214             else:
215                 self._mr_alt, self._mrName_alt,  self._mrInd_alt = r, Type, ind
216         if retAll:
217             return asscalar(copy(self._mr_alt)), self._mrName_alt, asscalar(copy(self._mrInd_alt))
218         else: return asscalar(copy(self._mr_alt))
219
220
221     def dmr(self, retAll = False):
222         # returns direction for max residual decrease
223         #( gradient for equality < 0 residuals ! )
224         return self.__dmr(retAll)
225
226     def __dmr(self, retAll = False):
227         if not hasattr(self, '_dmr') or (retAll and not hasattr(self, '_dmrInd')):
228             g = zeros(self.p.n)
229             maxResidual, resType, ind = self.mr(retAll=True)
230             if resType == 'lb':
231                 g[ind] -= 1 # N * (-1), -1 = dConstr/dx = d(lb-x)/dx
232             elif resType == 'ub':
233                 g[ind] += 1 # N * (+1), +1 = dConstr/dx = d(x-ub)/dx
234             elif resType == 'lin_ineq':
235                 g += self.p.A[ind]
236             elif resType == 'lin_eq':
237                 rr = self.p.matmult(self.p.Aeq[ind], self.x)-self.p.beq[ind]
238                 if rr < 0:  g -= self.p.Aeq[ind]
239                 else:  g += self.p.Aeq[ind]
240             elif resType == 'c':
241                 dc = self.dc(ind=ind).flatten()
242                 g += dc
243             elif resType == 'h':
244                 dh = self.dh(ind=ind).flatten()
245                 if self.p.h(self.x, ind) < 0:  g -= dh#CHECKME!!
246                 else: g += dh#CHECKME!!
247             else:
248                 # TODO: error or debug warning
249                 pass
250                 #self.p.err('incorrect resType')
251
252             self._dmr, self._dmrName,  self._dmrInd = g, resType, ind
253         if retAll:
254             return copy(self._dmr),  self._dmrName,  copy(self._dmrInd)
255         else:
256             return copy(self._dmr)
257
258     def betterThan(self, *args, **kwargs):
259         """
260         usage: result = involvedPoint.better(pointToCompare)
261
262         returns True if the involvedPoint is better than pointToCompare
263         and False otherwise
264         (if NOT better, mb same fval and same residuals or residuals less than desired contol)
265         """
266         return self.__betterThan__(*args, **kwargs)
267
268     def __betterThan__(self, point2compare, altLinInEq = False):
269         if self.p.isUC:
270             return self.f() < point2compare.f()
271
272         if altLinInEq:
273             mr_field = 'mr_alt'
274         else:
275             mr_field = 'mr'
276         point2compareResidual = getattr(point2compare, mr_field)()
277
278         criticalResidualValue = max((self.p.contol, point2compareResidual))
279
280         if hasattr(self, '_'+mr_field):
281             if getattr(self, '_'+mr_field) > criticalResidualValue: return False
282         else:
283             #TODO: simplify it!
284             #for fn in Residuals: (...)
285             if altLinInEq:
286                 if self.__all_lin_ineq() > criticalResidualValue: return False
287             else:
288                 if any(self.lb() > criticalResidualValue): return False
289                 if any(self.ub() > criticalResidualValue): return False
290                 if any(self.lin_ineq() > criticalResidualValue): return False
291                 if any(abs(self.lin_eq()) > criticalResidualValue): return False
292             if any(abs(self.h()) > criticalResidualValue): return False
293             if any(self.c() > criticalResidualValue): return False
294
295         mr = getattr(self, mr_field)()
296
297         if not self.p.isNaNInConstraintsAllowed:
298             if point2compare.__nNaNs__()  > self.__nNaNs__(): return True
299             elif point2compare.__nNaNs__()  < self.__nNaNs__(): return False
300             # TODO: check me
301             if mr <= self.p.contol and point2compareResidual <= self.p.contol and self.__nNaNs__() != 0: return mr < point2compareResidual
302
303         if mr < point2compareResidual and self.p.contol < point2compareResidual: return True
304
305         point2compareF_is_NaN = isnan(point2compare.f())
306         selfF_is_NaN = isnan(self.f())
307
308         if not point2compareF_is_NaN: # f(point2compare) is not NaN
309             if not selfF_is_NaN: # f(newPoint) is not NaN
310                 return self.f() < point2compare.f()
311             else: # f(newPoint) is NaN
312                 return False
313         else: # f(point2compare) is NaN
314             if selfF_is_NaN: # f(newPoint) is NaN
315                 return mr < point2compareResidual
316             else: # f(newPoint) is not NaN
317                 return True
318
319     def isFeas(self, **kwargs):
320         return self.__isFeas__(**kwargs)
321
322     def __isFeas__(self, altLinInEq = False):
323         if not all(isfinite(self.f())): return False
324         contol = self.p.contol
325         if altLinInEq:
326             if hasattr(self, '_mr_alt'):
327                 if self._mr_alt > contol or (not self.p.isNaNInConstraintsAllowed and self.__nNaNs__() != 0): return False
328             else:
329                 #TODO: simplify it!
330                 #for fn in Residuals: (...)
331                 if self.all_lin_ineq() > contol: return False
332         else:
333             if hasattr(self, '_mr'):
334                 if self._mr > contol or (not self.p.isNaNInConstraintsAllowed and self.__nNaNs__() != 0): return False
335             else:
336                 #TODO: simplify it!
337                 #for fn in Residuals: (...)
338                 if any(self.lb() > contol): return False
339                 if any(self.ub() > contol): return False
340                 if any(abs(self.lin_eq()) > contol): return False
341                 if any(self.lin_ineq() > contol): return False
342         if any(abs(self.h()) > contol): return False
343         if any(self.c() > contol): return False
344         return True
345
346     def __nNaNs__(self):
347         # returns number of nans in constraints
348         r = 0
349         c, h = self.c(), self.h()
350         r += len(where(isnan(c))[0])
351         r += len(where(isnan(h))[0])
352         return r
353
354     def directionType(self, *args, **kwargs):
355         return self.__directionType__(*args, **kwargs)
356
357     def __directionType__(self, *args, **kwargs):
358         if not hasattr(self, 'dType'):
359             self.__getDirection__(*args, **kwargs)
360         return self.dType
361
362     #def __getDirection__(self, useCurrentBestFeasiblePoint = False):
363     def __getDirection__(self, altLinInEq = False):
364         if hasattr(self, 'direction'):
365             return self.direction
366         p = self.p
367         contol = p.contol
368         maxRes, fname, ind = self.mr_alt(1)
369         if self.isFeas(altLinInEq=altLinInEq):
370         #or (useCurrentBestFeasiblePoint and hasattr(p, 'currentBestFeasiblePoint') and self.f() - p.currentBestFeasiblePoint.f() > self.mr()):
371         #if (maxRes <= p.contol and all(isfinite(self.df())) and (p.isNaNInConstraintsAllowed or self.__nNaNs__() == 0)) :
372             self.direction, self.dType = self.df(),'f'
373             return self.direction
374         else:
375             #d = zeros(p.n)
376             #if any(self.lb()>contol) or any(self.ub()>contol) or any(self.lin_eq()>contol) or any(self.lin_ineq()>contol):
377 #            lin_eq = self.lin_eq()
378 #            c = self.c()
379 #            h = self.h()
380
381 #            LB = lb[lb>contol/2]
382 #            UB = ub[ub>contol/2]
383 #            LIN_INEQ = lin_ineq[lin_ineq>contol/2]
384 #            nActiveLinInEq = LB.size + UB.size + LIN_INEQ.size
385 #            LinConstraints = sum(LB** 2) + sum(UB ** 2)
386 #            if  LIN_INEQ.size > 0: LinConstraints+= sum(LIN_INEQ ** 2)
387 #            #if  lin_eq.size > 0: LinConstraints+= sum(lin_eq[abs(lin_eq)>contol] ** 2)
388 ##
389 #            maxNonLinConstraint = 0.0
390 #            if c.size > 0: maxNonLinConstraint = max((maxNonLinConstraint, max(c)))
391 #            if h.size > 0: maxNonLinConstraint = max((maxNonLinConstraint, max(abs(h))))
392
393             if fname == 'all_lin_ineq':
394                 d = self.__all_lin_ineq_gradient()
395                 self.direction, self.dType = d, 'all_lin_ineq'
396             elif fname == 'lin_eq':
397                 raise "OpenOpt kernal error"
398                 #d = self.dmr()
399                 #self.dType = 'lin_eq'
400             elif fname == 'c':
401                 d = self.dmr()
402                 if p.debug: assert array_equal(self.dc(ind).flatten(), self.dmr())
403                 self.dType = 'c'
404             elif fname == 'h':
405                 d = self.dmr()#sign(self.h(ind))*self.dh(ind)
406                 if p.debug: assert array_equal(self.dh(ind).flatten(), self.dmr())
407                 self.dType = 'h'
408             else:
409                 p.err('error in getRalgDirection (unknown residual type ' + fname + ' ), you should report the bug')
410             self.direction = d.flatten()
411             return self.direction
412
413 #    def __getDirection__(self, useCurrentBestFeasiblePoint = False):
414 #        if hasattr(self, 'direction'):
415 #            return self.direction
416 #        p = self.p
417 #        contol = p.contol
418 #        maxRes, fname, ind = self.mr(retAll=True)
419 #        if (maxRes <= p.contol and all(isfinite(self.df())) and (p.isNaNInConstraintsAllowed or self.__nNaNs__() == 0)) \
420 #        or (useCurrentBestFeasiblePoint and hasattr(p, 'currentBestFeasiblePoint') and self.f() - p.currentBestFeasiblePoint.f() > self.mr()):
421 #        #if (maxRes <= p.contol and all(isfinite(self.df())) and (p.isNaNInConstraintsAllowed or self.__nNaNs__() == 0)) :
422 #            self.direction, self.dType = self.df(),'f'
423 #            return self.direction
424 #        else:
425 #            d = zeros(p.n)
426 #            #if any(self.lb()>contol) or any(self.ub()>contol) or any(self.lin_eq()>contol) or any(self.lin_ineq()>contol):
427 #            lb = self.lb()
428 #            ub = self.ub()
429 #            lin_ineq = self.lin_ineq()
430 #            lin_eq = self.lin_eq()
431 #            c = self.c()
432 #            h = self.h()
433 #
434 #            LB = lb[lb>contol/2]
435 #            UB = ub[ub>contol/2]
436 #            LIN_INEQ = lin_ineq[lin_ineq>contol/2]
437 #            nActiveLinInEq = LB.size + UB.size + LIN_INEQ.size
438 #            LinConstraints = sum(LB** 2) + sum(UB ** 2)
439 #            if  LIN_INEQ.size > 0: LinConstraints+= sum(LIN_INEQ ** 2)
440 #            #if  lin_eq.size > 0: LinConstraints+= sum(lin_eq[abs(lin_eq)>contol] ** 2)
441 ##
442 #            maxNonLinConstraint = 0.0
443 #            if c.size > 0: maxNonLinConstraint = max((maxNonLinConstraint, max(c)))
444 #            if h.size > 0: maxNonLinConstraint = max((maxNonLinConstraint, max(abs(h))))
445 #
446 #            if nActiveLinInEq != 0 and LinConstraints/nActiveLinInEq >= p.contol * maxNonLinConstraint and (lin_eq.size == 0 or LinConstraints/nActiveLinInEq >= p.contol * abs(lin_eq).max()):#fname in ['lb',  'ub',  'lin_eq',  'lin_ineq']:# or tmp > maxNonLinConstraint:
447 #                threshold = contol
448 #                ind_lb = where(lb>contol)[0]
449 #                ind_ub = where(ub>contol)[0]
450 #                ind_lin_ineq = where(lin_ineq>threshold)[0]
451 #                #ind_lin_eq = where(abs(lin_eq)>threshold)[0]
452 #
453 #                if ind_lb.size != 0:
454 #                    d[ind_lb] -= lb[ind_lb]# 0.5*d/dx((x-lb)^2) for violated constraints
455 #                if ind_ub.size != 0:
456 #                    d[ind_ub] += ub[ind_ub]# 0.5*d/dx((x-ub)^2) for violated constraints
457 #                if ind_lin_ineq.size != 0:
458 #                    a = p.A[ind_lin_ineq]
459 #                    b = p.b[ind_lin_ineq]
460 #                    d += dot(a.T, dot(a, self.x)  - b) # 0.5*d/dx((Ax-b)^2)
461 ##                if ind_lin_eq.size != 0:
462 ##                    aeq = p.Aeq[ind_lin_eq]
463 ##                    beq = p.beq[ind_lin_eq]
464 ##                    d += dot(aeq.T, dot(aeq, self.x)  - beq) # 0.5*d/dx((Ax-b)^2)
465 #                self.direction, self.dType = d/p.contol/nActiveLinInEq, 'linear'
466 #            elif fname == 'lin_eq':
467 #                d = self.dmr()
468 #                self.direction = d.flatten()
469 #                self.dType = 'lin_eq'
470 #            elif fname == 'c':
471 #                d = self.dmr()#self.dc(ind)
472 #                self.direction = d.flatten()
473 #                self.dType = 'c'
474 #            elif fname == 'h':
475 #                d = self.dmr()#sign(self.h(ind))*self.dh(ind)
476 #                self.direction = d.flatten()
477 #                self.dType = 'h'
478 #            else:
479 #                p.err('error in getRalgDirection (unknown residual type ' + fname + ' ), you should report the bug')
480 #            return self.direction
Note: See TracBrowser for help on using the browser.