Source Code for Module minfx.hessian_mods.gmw81_old

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003-2013 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the minfx optimisation library,                        # 
  6  # https://sourceforge.net/projects/minfx                                      # 
  7  #                                                                             # 
  8  # This program is free software: you can redistribute it and/or modify        # 
  9  # it under the terms of the GNU General Public License as published by        # 
 10  # the Free Software Foundation, either version 3 of the License, or           # 
 11  # (at your option) any later version.                                         # 
 12  #                                                                             # 
 13  # This program is distributed in the hope that it will be useful,             # 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 16  # GNU General Public License for more details.                                # 
 17  #                                                                             # 
 18  # You should have received a copy of the GNU General Public License           # 
 19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Module docstring. 
 24  """The Gill, Murray, and Wright (GMW) modified Cholesky Hessian modification algorithm. 
 25   
 26  Warning:  This code is not currently functional. 
 27   
 28  This file is part of the U{minfx optimisation library<https://sourceforge.net/projects/minfx>}. 
 29  """ 
 30   
 31  # Python module imports. 
 32  from numpy import dot, sort, sqrt, transpose 
 33  from numpy.linalg import eig, inv 
 34   
 35   
 36 -def gmw_old(dfk, d2fk, I, n, mach_acc, print_prefix, print_flag, return_matrix=0): 
 37      """Modified Cholesky Hessian modification. 
 38   
 39      Algorithm 6.5 from page 148. 
 40   
 41      Somehow the data structures l, d, and c can be stored in d2fk! 
 42      """ 
 43   
 44      # Debugging (REMOVE!!!). 
 45      #d2fk = array([[4, 2, 1], [2, 6, 3], [1, 3, -0.004]], float64) 
 46   
 47      if print_flag >= 3: 
 48          print("d2fk: " + repr(d2fk)) 
 49          eigen = eig(d2fk) 
 50          eigenvals = sort(eigen[0]) 
 51          print("Eigenvalues: " + repr(eigenvals)) 
 52   
 53      # Calculate gamma(A) and xi(A). 
 54      gamma = 0.0 
 55      xi = 0.0 
 56      for i in range(n): 
 57          gamma = max(abs(d2fk[i, i]), gamma) 
 58          for j in range(i+1, n): 
 59              xi = max(abs(d2fk[i, j]), xi) 
 60   
 61      # Calculate delta and beta. 
 62      delta = mach_acc * max(gamma + xi, 1) 
 63      beta = sqrt(max(gamma, xi / sqrt(n**2 - 1.0), mach_acc)) 
 64   
 65      # Initialise data structures. 
 66      d2fk_orig = 1.0 * d2fk 
 67      c = 0.0 * d2fk 
 68      d = 0.0 * dfk 
 69      l = 1.0 * I 
 70      e = 0.0 * dfk 
 71      P = 1.0 * I 
 72   
 73      # Initial diagonal elements of c. 
 74      for k in range(n): 
 75          c[k, k] = d2fk[k, k] 
 76   
 77      # Main loop. 
 78      for j in range(n): 
 79          if print_flag >= 3: 
 80              print("\n<j: " + repr(j) + ">") 
 81   
 82          # Row and column swapping. 
 83          p = 1.0 * I 
 84          q = j 
 85          for i in range(j, n): 
 86              if abs(c[q, q]) <= abs(c[i, i]): 
 87                  q = i 
 88          if print_flag >= 3: 
 89              print("Row and column swapping.") 
 90              print("i range: " + repr(range(j, n))) 
 91              print("q: " + repr(q)) 
 92              print("a: " + repr(d2fk)) 
 93              print("c: " + repr(c)) 
 94              print("d: " + repr(d)) 
 95              print("l: " + repr(l)) 
 96          if q != j: 
 97              # Modify the permutation matrices. 
 98              temp_p, temp_P = 1.0*p[:, q], 1.0*P[:, q] 
 99              p[:, q], P[:, q] = p[:, j], P[:, j] 
100              p[:, j], P[:, j] = temp_p, temp_P 
101   
102              # Permute both c and d2fk. 
103              c = dot(p, dot(c, p)) 
104              d2fk = dot(p, dot(d2fk, p)) 
105   
106              if print_flag >= 3: 
107                  print("p: " + repr(p)) 
108                  print("a(mod): " + repr(d2fk)) 
109                  print("c(mod): " + repr(c)) 
110                  print("d(mod): " + repr(d)) 
111                  print("l(mod): " + repr(l)) 
112   
113          # Calculate the elements of l. 
114          if print_flag >= 3: 
115              print("\nCalculate the elements of l.") 
116              print("s range: " + repr(range(j))) 
117          for s in range(j): 
118              if print_flag >= 3: 
119                  print("s: " + repr(s)) 
120              l[j, s] = c[j, s] / d[s] 
121          if print_flag >= 3: 
122              print("l: " + repr(l)) 
123   
124          # Calculate c[i, j]. 
125          if print_flag >= 3: 
126              print("\nCalculate c[i, j].") 
127              print("i range: " + repr(range(j+1, n))) 
128          for i in range(j+1, n): 
129              sum = 0.0 
130              for s in range(j): 
131                  if print_flag >= 3: 
132                      print("s range: " + repr(range(j))) 
133                  sum = sum + l[j, s] * c[i, s] 
134              c[i, j] = d2fk[i, j] - sum 
135          if print_flag >= 3: 
136              print("c: " + repr(c)) 
137   
138          # Calculate d. 
139          if print_flag >= 3: 
140              print("\nCalculate d.") 
141          theta_j = 0.0 
142          if j < n-1: 
143              if print_flag >= 3: 
144                  print("j < n, " + repr(j) + " < " + repr(n-1)) 
145                  print("i range: " + repr(range(j+1, n))) 
146              for i in range(j+1, n): 
147                  theta_j = max(theta_j, abs(c[i, j])) 
148          else: 
149              if print_flag >= 3: 
150                  print("j >= n, " + repr(j) + " >= " + repr(n-1)) 
151          d[j] = max(abs(c[j, j]), (theta_j/beta)**2, delta) 
152          if print_flag >= 3: 
153              print("theta_j: " + repr(theta_j)) 
154              print("d: " + repr(d)) 
155   
156          # Calculate c[i, i]. 
157          if print_flag >= 3: 
158              print("\nCalculate c[i, i].") 
159          if j < n-1: 
160              if print_flag >= 3: 
161                  print("j < n, " + repr(j) + " < " + repr(n-1)) 
162                  print("i range: " + repr(range(j+1, n))) 
163              for i in range(j+1, n): 
164                  c[i, i] = c[i, i] - c[i, j]**2 / d[j] 
165          else: 
166              if print_flag >= 3: 
167                  print("j >= n, " + repr(j) + " >= " + repr(n-1)) 
168          if print_flag >= 3: 
169              print("c: " + repr(c)) 
170   
171          # Calculate e. 
172          e[j] = d[j] - c[j, j] 
173   
174      for i in range(n): 
175          d2fk[i, i] = d2fk[i, i] + e[i] 
176      d2fk = dot(P, dot(d2fk, transpose(P))) 
177   
178      # Debugging. 
179      if print_flag >= 3: 
180          print("\nFin:") 
181          print("gamma(A): " + repr(gamma)) 
182          print("xi(A): " + repr(xi)) 
183          print("delta: " + repr(delta)) 
184          print("beta: " + repr(beta)) 
185          print("c: " + repr(c)) 
186          print("d: " + repr(d)) 
187          print("l: " + repr(l)) 
188          temp = 0.0 * I 
189          for i in range(len(d)): 
190              temp[i, i] = d[i] 
191          print("m: " + repr(dot(l, sqrt(temp)))) 
192          print("mmT: " + repr(dot(dot(l, sqrt(temp)), transpose(dot(l, sqrt(temp)))))) 
193          print("P: " + repr(P)) 
194          print("e: " + repr(e)) 
195          print("d2fk: " + repr(d2fk)) 
196          eigen = eig(d2fk) 
197          eigenvals = sort(eigen[0]) 
198          print("Eigenvalues: " + repr(eigenvals)) 
199          import sys 
200          sys.exit() 
201   
202      # Calculate the Newton direction. 
203      if return_matrix: 
204          return -dot(inv(d2fk), dfk), d2fk 
205      else: 
206          return -dot(inv(d2fk), dfk) 
207