Source Code for Module parasol.Goal

  1   
  2   
  3  __author__ = "Charlie Taylor (charlietaylor@sourceforge.net)" 
  4  __version__ = " 1.0 " 
  5  __date__ = "Apr 2, 2005" 
  6  __copyright__ = "Copyright (c) 2009 Charlie Taylor" 
  7  __license__ = "BSD" 
  8   
  9  from math import * 
 10   
 11 -class Goal( object ): 
 12       
 13      """ 
 14       Find the value of X in the range minX to maxX that returns goalVal from funcOfX. 
 15       
 16       This CLASS version of goal finder allows nested Goal loops, 
 17       (unlike FORTRAN version via f2py) 
 18   
 19       Based on the world's best root finder by Jack Crenshaw, 
 20       U{http://www.embedded.com/story/OEG20030508S0030}, 
 21       turned into object-oriented world's best "goal-value" finder by Charlie Taylor 
 22        
 23       returns (X, ierror): the value of X and an error flag 
 24   
 25      """ 
 26           
 27 -    def __call__(self, xval=0.0): 
 28           
 29          ierror = 1 # needs to be reset in try block to avoid fail logic 
 30          try: 
 31              X, ierror = self.cren_goal(self.goalVal, self.minX, self.maxX,  
 32                  self.funcOfX, self.tolerance, self.maxLoops) 
 33          except: 
 34              print "ERROR unrecoverable in Goal.py" 
 35              X = self.maxX 
 36               
 37          if ierror: 
 38              if (self.failValue==None): 
 39                  X = self.maxX 
 40                  print "returned an ERROR=%s from cren_goal in goalFind.py"%str(ierror) 
 41                  print "goalVal, minX, maxX, funcOfX",self.goalVal, self.minX, self.maxX, self.funcOfX.__name__ 
 42                  print "using X value=",X 
 43              else: 
 44                  X = self.failValue 
 45               
 46          return X, ierror 
 47           
 48   
 49          
 50 -    def __init__(self, goalVal=0.0, minX=0.0, maxX=1000.0,  
 51          funcOfX=None, tolerance=1.0E-4, maxLoops=40, failValue=None): 
 52               
 53          '''Find the value of X in the range minX to maxX that returns goalVal from funcOfX. 
 54              
 55             @param goalVal: goal value (float) 
 56             @param minX: minimum value of X considered (float) 
 57             @param maxX: maximum value of X considered (float) 
 58             @param funcOfX: function of x (callable),  
 59                 expects a float input, returns a float.  
 60             @param tolerance: accuracy termination criteria (float), in same units as goalVal 
 61             @param maxLoops: maximum number of search loops (int) 
 62             @param failValue: value returned in case of error (float or None) 
 63          ''' 
 64   
 65        # initialize properties 
 66          self.goalVal = goalVal 
 67          self.minX = minX 
 68          self.maxX = maxX 
 69          self.funcOfX = funcOfX 
 70          self.tolerance = tolerance 
 71          self.maxLoops = maxLoops 
 72          self.failValue = failValue 
 73       
 74       
 75 -    def cren_goal(self, goal, x0, x2, f, eps, imax): 
 76          # ierror parameters 
 77          SUCCESS = 0 
 78          MyBAD_DATA = 1 
 79          NO_CONVERGENCE = 2 
 80             
 81          xmlast = x0 
 82           
 83          y0 = f(x0) - goal 
 84          if (y0 == 0.0): 
 85              return x0, SUCCESS 
 86         
 87          y2 = f(x2) - goal 
 88          if (y2 == 0.0): 
 89              return x2, SUCCESS 
 90       
 91          if(y2 * y0 > 0.0): 
 92              return x0, MyBAD_DATA 
 93           
 94          for i in range(imax): 
 95              #print'I=',i,' of ',imax 
 96             
 97              x1 = 0.5 * (x2 + x0) 
 98              y1 = f(x1) - goal 
 99              if (y1 ==  0.0) : 
100                  return x1, SUCCESS 
101               
102              if (abs (x1 - x0) < eps) : 
103                  return x1, SUCCESS 
104               
105              if (y1 * y0 > 0.0): 
106                  x0,x2,y0,y2 = (x2,x0,y2,y0) # switch points 0 and 2 
107               
108              y10 = y1 - y0 
109              y21 = y2 - y1 
110              y20 = y2 - y0 
111              if (y2 * y20 < 2.0 * y1 * y10): 
112                  x2,y2 = (x1,y1) 
113              else: 
114                  b = (x1 - x0) / y10    
115                  c = (y10 -y21) / (y21 * y20)  
116                  xm = x0 - b * y0 * (1.0-c * y1) 
117                  ym = f(xm) - goal 
118                  if (ym ==  0.0)  : 
119                      return xm, SUCCESS 
120                   
121                  if (abs (xm - xmlast) < eps) : 
122                      return xm, SUCCESS 
123                   
124                  xmlast = xm 
125                  if (ym * y0 < 0.0): 
126                      x2,y2 = (xm,ym) 
127                  else: 
128                      x0,y0,x2,y2 = (xm,ym,x1,y1) 
129                   
130          return x1, NO_CONVERGENCE 
131               
132   
133  if __name__ == "__main__":  #self test 
134       
135 -    def FofX(x): 
136          #print 'calling FofX with x=',x 
137          return x**3 - x**2 + x - 11.875 
138           
139      print '======================================================' 
140      print '----->Answer should be = 2.5' 
141      print 'with tolerance=1.0E-3  =', 
142      G = Goal(goalVal=0.0, minX=0.0, maxX=57.123,  
143          funcOfX=FofX, tolerance=1.0E-3, maxLoops=40, failValue=None) 
144      print G() 
145       
146      print 'with tolerance=1.0E-20 =', 
147      G2 = Goal(goalVal=0.0, minX=0.0, maxX=57.123,  
148          funcOfX=FofX, tolerance=1.0E-20, maxLoops=40, failValue=None) 
149      print G2() 
150