Changeset 828 for mystic

Timestamp:

09/18/15 16:05:53 (8 months ago)

Author:

mmckerns

Message:

update SV regression examples from examples_other

Location:

mystic

Files:

• examples/test_svr1.py (moved) (moved from mystic/examples_other/test_svr1.py) (5 diffs)
• examples/test_svr2.py (moved) (moved from mystic/examples_other/test_svr2.py) (3 diffs)
• examples_other/qld_svr1.py (copied) (copied from mystic/examples_other/test_svr1.py)
• examples_other/qld_svr2.py (copied) (copied from mystic/examples_other/test_svr2.py)

mystic/examples/test_svr1.py

@@ -2 +2 @@
 #
 # Author: Patrick Hung (patrickh @caltech)
+# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
 # Copyright (c) 1997-2015 California Institute of Technology.
 # License: 3-clause BSD.  The full license text is available at:
@@ -10 +11 @@
 
 from numpy import *
-import qld
 import pylab
 from mystic.svmtools import *
 
-from mystic.tools import random_seed
-#random_seed(123)
+# a common objective function for solving a QP problem
+# (see http://www.mathworks.com/help/optim/ug/quadprog.html)
+def objective(x, H, f):
+    return 0.5 * dot(dot(x,H),x) + dot(f,x)
 
+# linear data with scatter
 x = arange(-5, 5.001)
-y = x + random.rand(x.size)*7
-#y = array([9.99, 10.93, 9.11, 6.50, 12.69, 14.09, 13.83, 15.86, 17.49, 14.20, 19.98])
+y = x + 7*random.rand(x.size)
 
 l = x.size
 N = 2*l
-
-pylab.plot(x, y, 'k+')
-#pylab.show()
 
 # the Kernel Matrix (with the linear kernel)
@@ -33 +32 @@
 # the linear term for the QP
 Y = concatenate([y,-y])
-svr_epsilon = 3;
+svr_epsilon = 3
 b = Y + svr_epsilon * ones(Y.size)
 
@@ -39 +38 @@
 f = b
 Aeq = concatenate([ones(l), -ones(l)]).reshape(1,N)
-Beq = array([0])
+Beq = array([0.])
 lb = zeros(N)
 ub = zeros(N) + 0.5
 
+from mystic.symbolic import linear_symbolic, solve, \
+     generate_solvers as solvers, generate_constraint as constraint
+constrain = linear_symbolic(Aeq,Beq)
+constrain = constraint(solvers(solve(constrain,target=['x0'])))
 
-alpha = qld.quadprog2(H, f, None, None, Aeq, Beq, lb, ub)
-#alpha = 2*alpha
+from mystic import supressed
+@supressed(1e-5)
+def conserve(x):
+    return constrain(x)
+
+from mystic.monitors import VerboseMonitor
+mon = VerboseMonitor(10)
+
+from mystic.solvers import diffev
+alpha = diffev(objective, zip(lb,.1*ub), args=(H,f), npop=N*3, gtol=200, \
+               itermon=mon, \
+               ftol=1e-5, bounds=zip(lb,ub), constraints=conserve, disp=1)
+
+print 'solved x: ', alpha
+print "constraint A*x == 0: ", inner(Aeq, alpha)
+print "minimum 0.5*x'Hx + f'*x: ", objective(alpha, H, f)
 
 sv1 = SupportVectors(alpha[:l])
@@ -52 +69 @@
 R = RegressionFunction(x, y, alpha, svr_epsilon, LinearKernel)
 
-pylab.plot(x,R(x))
-pylab.plot(x,R(x)+svr_epsilon, 'k--')
-pylab.plot(x,R(x)-svr_epsilon, 'k--')
+print 'support vectors: ', sv1, sv2
+
+# plot data
+pylab.plot(x, y, 'k+', markersize=10)
+
+# plot regression function and support
+pylab.plot(x,R(x), 'k-')
+pylab.plot(x,R(x)-svr_epsilon, 'r--')
+pylab.plot(x,R(x)+svr_epsilon, 'g--')
 pylab.plot(x[sv1],y[sv1],'ro')
 pylab.plot(x[sv2],y[sv2],'go')
-
-print alpha
-
 pylab.show()
 
-# $Id$
-# 
 # end of file

mystic/examples/test_svr2.py

@@ -2 +2 @@
 #
 # Author: Patrick Hung (patrickh @caltech)
+# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
 # Copyright (c) 1997-2015 California Institute of Technology.
 # License: 3-clause BSD.  The full license text is available at:
 #  - http://mmckerns.github.io/project/mystic/browser/mystic/LICENSE
 """
-Support Vector Regression. Example 1
+Support Vector Regression. Example 2
 """
 
 from numpy import *
-import qld
 import pylab
 from mystic.svmtools import *
 
+# a common objective function for solving a QP problem
+# (see http://www.mathworks.com/help/optim/ug/quadprog.html)
+def objective(x, H, f):
+    return 0.5 * dot(dot(x,H),x) + dot(f,x)
+
+# quadratic data plus scatter
 x = arange(-5, 5.001)
-y = x*x + random.rand(x.size)*8
-#y = array([9.99, 10.93, 9.11, 6.50, 12.69, 14.09, 13.83, 15.86, 17.49, 14.20, 19.98])
+y = x*x + 8*random.rand(x.size)
 
 l = x.size
 N = 2*l
-
-pylab.plot(x, y, 'k+')
 
 # the Kernel Matrix 
@@ -37 +40 @@
 f = b
 Aeq = concatenate([ones(l), -ones(l)]).reshape(1,N)
-Beq = array([0])
+Beq = array([0.])
 lb = zeros(N)
 ub = zeros(N) + 0.5
 
-alpha = qld.quadprog2(H, f, None, None, Aeq, Beq, lb, ub)
-#alpha = 2*alpha
+from mystic.symbolic import linear_symbolic, solve, \
+     generate_solvers as solvers, generate_constraint as constraint
+constrain = linear_symbolic(Aeq,Beq)
+constrain = constraint(solvers(solve(constrain,target=['x0'])))
+
+from mystic import supressed
+@supressed(1e-5)
+def conserve(x):
+    return constrain(x)
+
+from mystic.monitors import VerboseMonitor
+mon = VerboseMonitor(10)
+
+from mystic.solvers import diffev
+alpha = diffev(objective, zip(lb,.1*ub), args=(H,f), npop=N*3, gtol=200, \
+               itermon=mon, \
+               ftol=1e-5, bounds=zip(lb,ub), constraints=conserve, disp=1)
+
+print 'solved x: ', alpha
+print "constraint A*x == 0: ", inner(Aeq, alpha)
+print "minimum 0.5*x'Hx + f'*x: ", objective(alpha, H, f)
 
 sv1 = SupportVectors(alpha[:l])
@@ -49 +71 @@
 R = RegressionFunction(x, y, alpha, svr_epsilon, pk)
 
+print 'support vectors: ', sv1, sv2
+
+# plot data
+pylab.plot(x, y, 'k+', markersize=10)
+
+# plot regression function and support
 xx = arange(min(x),max(x),0.1)
 pylab.plot(xx,R(xx))
-pylab.plot(xx,R(xx)+svr_epsilon, 'k--')
-pylab.plot(xx,R(xx)-svr_epsilon, 'k--')
+pylab.plot(xx,R(xx)-svr_epsilon, 'r--')
+pylab.plot(xx,R(xx)+svr_epsilon, 'g--')
 pylab.plot(x[sv1],y[sv1],'ro')
 pylab.plot(x[sv2],y[sv2],'go')
-
-print alpha
-
 pylab.show()
 
-# $Id$
-# 
 # end of file