Changeset 830 for mystic/examples/test_svc2.py

Timestamp:

09/19/15 11:26:58 (8 months ago)

Author:

mmckerns

Message:

cleaning and documentation for SVR and SVC examples; rename StandardLinearKernel?

File:

• mystic/examples/test_svc2.py (modified) (6 diffs)

mystic/examples/test_svc2.py

@@ -17 +17 @@
 import os.path
 
-# 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)
+# define the objective function as the dual for SVC
+# (see: http://scikit-learn.org/stable/modules/svm.html#svc)
+# the objective funciton is very similar to the dual for SVC
+# (see: http://scikit-learn.org/stable/modules/svm.html#svc)
+def objective(x, Q, b):
+    return 0.5 * dot(dot(x,Q),x) + dot(b,x)
 
 # SETTINGS
@@ -26 +28 @@
 overlap = False # reduce the distance between the datasets
 
+# define the data points for each class
 c1 = loadtxt(os.path.join('DATA','g1.pts'))
 c2 = loadtxt(os.path.join('DATA','g2.pts'))
@@ -34 +37 @@
 if overlap: c1[:,0] += 3 # make the datasets overlap a little
 
-# the Kernel Matrix (with the linear kernel)
-# Q = multiply.outer(X,X)   # NOTE: requires X is a list of scalars
+# define the full data set
+X = concatenate([c1,c2]); nx = X.shape[0]
+# define the labels (+1 for c1; -1 for c2)
+y = concatenate([ones(c1.shape[0]), -ones(c2.shape[0])]).reshape(1,nx)
 
+# build the Kernel Matrix (with the linear kernel)
+# get the QP quadratic and linear terms
 XX = concatenate([c1,-c2])
-nx = XX.shape[0]
+Q = KernelMatrix(XX)  # Q_ij = K(x_i, x_j)
+b = -ones(nx)         # b_i = 1  (in dual)
 
-# quadratic and linear terms of QP
-Q = KernelMatrix(XX)
-b = -1 * ones(nx)
-
-H = Q
-f = b
-Aeq = concatenate([ones(c1.shape[0]), -ones(c2.shape[0])]).reshape(1,nx)
+# build the constraints (y.T * x = 0.0)
+# 1.0*x0 + 1.0*x1 + ... - 1.0*xN = 0.0
+Aeq = y
 Beq = array([0.])
+# set the bounds
 lb = zeros(nx)
 ub = 1 * ones(nx)
 _b = .1 * ones(nx) # good starting value if most solved xi should be zero
 
+# build the constraints operator
 from mystic.symbolic import linear_symbolic, solve, \
      generate_solvers as solvers, generate_constraint as constraint
@@ -71 +77 @@
 mon = VerboseMonitor(10)
 
+# solve the dual for alpha
 from mystic.solvers import diffev
-alpha = diffev(objective, zip(lb,_b), args=(H,f), npop=nx*3, gtol=200,\
+alpha = diffev(objective, zip(lb,_b), args=(Q,b), npop=nx*3, gtol=200,\
                itermon=mon, \
                ftol=1e-8, bounds=zip(lb,ub), constraints=conserve, disp=1)
@@ -78 +85 @@
 print 'solved x: ', alpha
 print "constraint A*x == 0: ", inner(Aeq, alpha)
-print "minimum 0.5*x'Hx + f'*x: ", objective(alpha, H, f)
+print "minimum 0.5*x'Qx + b'*x: ", objective(alpha, Q, b)
 
-# the labels and the points
-X = concatenate([c1,c2])
-y = concatenate([ones(c1.shape[0]), -ones(c2.shape[0])]).reshape(1,nx)
-
+# calculate weight vectors, support vectors, and bias
 wv = WeightVector(alpha, X, y)
 sv1, sv2 = SupportVectors(alpha,y,eps=1e-6)
@@ -110 +114 @@
 
 # plot the support points
-pylab.plot(XX[sv1,0], XX[sv1,1], 'b^',markersize=10)
-pylab.plot(-XX[sv2,0], -XX[sv2,1], 'y^',markersize=10)
+pylab.plot(XX[sv1,0], XX[sv1,1], 'bo', markersize=8)
+pylab.plot(-XX[sv2,0], -XX[sv2,1], 'yo', markersize=8)
 #pylab.axis('equal')
 pylab.show()