Source Code for Module minimise.line_search.interpolate

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003 Edward d'Auvergne                                        # 
  4  #                                                                             # 
  5  # This file is part of the program relax.                                     # 
  6  #                                                                             # 
  7  # relax is free software; you can redistribute it and/or modify               # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation; either version 2 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # relax is distributed in the hope that it will be useful,                    # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with relax; if not, write to the Free Software                        # 
 19  # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA   # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23   
 24  from Numeric import sqrt 
 25   
 26 -def cubic_int(a, b, fa, fb, ga, gb): 
 27      """Cubic interpolation using f(a), f(b), g(a), and g(b). 
 28   
 29      Equations 
 30      ~~~~~~~~~ 
 31   
 32      f(a) = a'a**3 + b'a**2 + c'a + d' 
 33      f(b) = a'b**3 + b'b**2 + c'b + d' 
 34      g(a) = 3a'a**2 + 2b'a + c' 
 35      g(b) = 3a'b**2 + 2b'b + c' 
 36   
 37   
 38      Interpolation 
 39      ~~~~~~~~~~~~~ 
 40   
 41      The extrema are the roots of the quadratic equation: 
 42   
 43          3a'*alpha**2 + 2b'*alpha + c' = 0 
 44   
 45      The cubic interpolant is given by the formula: 
 46   
 47                             g(b) + beta2 - beta1 
 48          ac = b - (b - a) . --------------------- 
 49                             g(b) - g(a) + 2*beta2 
 50   
 51      where: 
 52                                    f(a) - f(b) 
 53          beta1 = g(a) + g(b) - 3 . ----------- 
 54                                       a - b 
 55   
 56          if a < b: 
 57              beta2 = sqrt(beta1**2 - g(a).g(b)) 
 58          else: 
 59              beta2 = -sqrt(beta1**2 - g(a).g(b)) 
 60      """ 
 61   
 62   
 63      beta1 = ga + gb - 3.0*(fa - fb)/(a - b) 
 64      if a < b: 
 65          beta2 = sqrt(beta1**2 - ga*gb) 
 66      else: 
 67          beta2 = -sqrt(beta1**2 - ga*gb) 
 68   
 69      denom = gb - ga + 2.0*beta2 
 70      if denom == 0.0: 
 71          return -1e99 
 72      else: 
 73          return b - (b - a)*(gb + beta2 - beta1)/denom 
 74   
 75   
 76 -def cubic_ext(a, b, fa, fb, ga, gb, full_output=0): 
 77      """Cubic Extrapolation using f(a), f(b), g(a), and g(b). 
 78   
 79      Extrapolation 
 80      ~~~~~~~~~~~~~ 
 81   
 82      The extrema are the roots of the quadratic equation: 
 83   
 84          3a'*alpha**2 + 2b'*alpha + c' = 0 
 85   
 86      The cubic extrapolant is given by the formula: 
 87   
 88                             g(b) + beta2 - beta1 
 89          ac = b - (b - a) . --------------------- 
 90                             g(b) - g(a) + 2*beta2 
 91   
 92      where: 
 93   
 94                                    f(a) - f(b) 
 95          beta1 = g(a) + g(b) - 3 . ----------- 
 96                                       a - b 
 97   
 98          if a < b: 
 99              beta2 = sqrt(max(0.0, beta1**2 - g(a).g(b))) 
100          else: 
101              beta2 = -sqrt(max(0.0, beta1**2 - g(a).g(b))) 
102      """ 
103   
104   
105      beta1 = ga + gb - 3.0*(fa - fb)/(a - b) 
106      if a < b: 
107          beta2 = sqrt(max(0.0, beta1**2 - ga*gb)) 
108      else: 
109          beta2 = -sqrt(max(0.0, beta1**2 - ga*gb)) 
110   
111      denom = gb - ga + 2.0*beta2 
112      if denom == 0.0: 
113          alpha = -1e99 
114      else: 
115          alpha = b - (b - a)*(gb + beta2 - beta1)/denom 
116   
117      if full_output: 
118          return alpha, beta1, beta2 
119      else: 
120          return alpha 
121   
122   
123 -def quadratic_fafbga(a, b, fa, fb, ga): 
124      """Quadratic interpolation using f(a), f(b), and g(a). 
125   
126      The extremum of the quadratic is given by: 
127   
128                       1             g(a) 
129          aq  =  a  +  - . ------------------------- 
130                       2   f(a) - f(b) - (a - b)g(a) 
131      """ 
132   
133      denom = fa - fb - (a - b)*ga 
134      if denom == 0.0: 
135          return 1e99 
136      else: 
137          return a + 0.5 * ga / denom 
138   
139   
140 -def quadratic_gagb(a, b, ga, gb): 
141      """Quadratic interpolation using g(a) and g(b). 
142   
143      The extremum of the quadratic is given by: 
144   
145                 bg(a) - ag(b) 
146          aq  =  ------------- 
147                  g(a) - g(b) 
148      """ 
149   
150      denom = ga - gb 
151      if denom == 0.0: 
152          return 1e99 
153      else: 
154          return (b*ga - a*gb) / denom 
155