From 734cdcebfe66ee283e7ae535cec13a857f8b4a29 Mon Sep 17 00:00:00 2001
From: "trevor.cohn" <trevor.cohn@ec762483-ff6d-05da-a07a-a48fb63a330f>
Date: Sat, 10 Jul 2010 16:49:23 +0000
Subject: Testing code for line search with simplex *equality* constraints and
 positivity constraints.

git-svn-id: https://ws10smt.googlecode.com/svn/trunk@217 ec762483-ff6d-05da-a07a-a48fb63a330f
---
 gi/posterior-regularisation/simplex_pg.py | 55 +++++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)
 create mode 100644 gi/posterior-regularisation/simplex_pg.py

(limited to 'gi/posterior-regularisation/simplex_pg.py')

diff --git a/gi/posterior-regularisation/simplex_pg.py b/gi/posterior-regularisation/simplex_pg.py
new file mode 100644
index 00000000..5da796d3
--- /dev/null
+++ b/gi/posterior-regularisation/simplex_pg.py
@@ -0,0 +1,55 @@
+#
+# Following Leunberger and Ye, Linear and Nonlinear Progamming, 3rd ed. p367
+# "The gradient projection method"
+# applied to an equality constraint for a simplex.
+#
+#   min f(x)
+#   s.t. x >= 0, sum_i x = d
+#
+# FIXME: enforce the positivity constraint - a limit on the line search?
+#
+
+from numpy import *
+from scipy import *
+from linesearch import line_search
+# local copy of scipy's Amijo line_search - wasn't enforcing alpha max correctly
+import sys
+
+dims = 4
+
+def f(x):
+    fv = x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[3]*x[3] - 2*x[0] - 3*x[3]
+    # print 'evaluating f at', x, 'value', fv
+    return fv
+
+def g(x):
+    return array([2*x[0] - 2, 2*x[1], 2*x[2], 2*x[3]-3])
+
+def pg(x):
+    gv = g(x)
+    return gv - sum(gv) / dims
+
+x = ones(dims) / dims
+old_fval = None
+
+while True:
+    fv = f(x)
+    gv = g(x)
+    dv = pg(x)
+
+    print 'x', x, 'f', fv, 'g', gv, 'd', dv
+
+    if old_fval == None:
+        old_fval = fv + 0.1
+
+    # solve for maximum step size i.e. when positivity constraints kick in
+    # x - alpha d = 0   => alpha = x/d
+    amax = max(x/dv)
+    if amax < 1e-8: break
+
+    stuff = line_search(f, pg, x, -dv, dv, fv, old_fval, amax=amax)
+    alpha = stuff[0] # Nb. can avoid next evaluation of f,g,d using 'stuff'
+    if alpha < 1e-8: break
+    x -= alpha * dv
+
+    old_fval = fv
-- 
cgit v1.2.3