Changeset 827 for mystic


Ignore:
Timestamp:
09/18/15 09:51:13 (8 months ago)
Author:
mmckerns
Message:

fix SV classification examples, move from examples_other to examples

Location:
mystic/examples
Files:
2 copied
1 moved

Legend:

Unmodified
Added
Removed

  • mystic/examples/test_svc1.py

    r826 r827  
    22# 
    33# Author: Patrick Hung (patrickh @caltech) 
     4# Author: Mike McKerns (mmckerns @caltech and @uqfoundation) 
    45# Copyright (c) 1997-2015 California Institute of Technology. 
    56# License: 3-clause BSD.  The full license text is available at: 
    … …  
    1011 
    1112from numpy import * 
    12 import qld, quapro 
    1313import pylab 
    1414from mystic.svctools import * 
    1515 
    16 def myshow(): 
    17     import Image, tempfile 
    18     name = tempfile.mktemp() 
    19     pylab.savefig(name,dpi=72) 
    20     im = Image.open('%s.png' % name) 
    21     im.show() 
     16# a common objective function for solving a QP problem 
     17# (see http://www.mathworks.com/help/optim/ug/quadprog.html) 
     18def objective(x, H, f): 
     19    return 0.5 * dot(dot(x,H),x) + dot(f,x) 
    2220 
    2321c1 = array([[0., 0.],[1., 0.],[ 0.2, 0.2],[0.,1.]]) 
    2422c2 = array([[0, 1.1], [1.1, 0.],[0, 1.5],[0.5,1.2],[0.8, 1.7]]) 
    2523 
    26 pylab.plot(c1[:,0], c1[:,1], 'ko') 
    27 pylab.plot(c2[:,0], c2[:,1], 'ro') 
    28  
    2924# the Kernel Matrix (with the linear kernel) 
    30 # Q = multiply.outer(X,X) <--- unfortunately only works when X is a list of scalars... 
    31 # In Mathematica, this would be implemented simply via Outer[K, X, X, 1] 
     25# Q = multiply.outer(X,X)   # NOTE: requires X is a list of scalars 
    3226 
    3327XX = concatenate([c1,-c2]) 
    … …  
    4135f = b 
    4236Aeq = concatenate([ones(c1.shape[0]), -ones(c2.shape[0])]).reshape(1,nx) 
    43 Beq = array([0]) 
     37Beq = array([0.]) 
    4438lb = zeros(nx) 
    45 ub = zeros(nx) + 99999 
     39ub = 99999 * ones(nx) 
    4640 
    47 #alpha = qld.quadprog2(H, f, None, None, Aeq, Beq, lb, ub) 
    48 alpha, xxx= quapro.quadprog(H, f, None, None, Aeq=Aeq, beq=Beq, LB=lb, UB=ub) 
    49 print alpha 
    50 #alpha = array([fil(x) for x in alpha]) 
    51 print "cons:", inner(Aeq,alpha) 
    52 print "obj min: 0.5 * x'Hx + <x,f>",  0.5*inner(alpha,inner(H,alpha))+inner(f,alpha) 
     41from mystic.symbolic import linear_symbolic, solve, \ 
     42     generate_solvers as solvers, generate_constraint as constraint 
     43constrain = linear_symbolic(Aeq,Beq) 
     44constrain = constraint(solvers(solve(constrain,target=['x0']))) 
    5345 
     46from mystic import supressed 
     47@supressed(1e-5) 
     48def conserve(x): 
     49    return constrain(x) 
     50 
     51#from mystic.monitors import VerboseMonitor 
     52#mon = VerboseMonitor(1) 
     53 
     54from mystic.solvers import diffev 
     55alpha = diffev(objective, zip(lb,ub), args=(H,f), npop=nx*3, gtol=200, \ 
     56#              itermon=mon, \ 
     57               ftol=1e-8, bounds=zip(lb,ub), constraints=conserve, disp=1) 
     58 
     59print 'solved x: ', alpha 
     60print "constraint A*x == 0: ", inner(Aeq, alpha) 
     61print "minimum 0.5*x'Hx + f'*x: ", objective(alpha, H, f) 
    5462 
    5563# the labels and the points 
    … …  
    6674bias2 = -0.5 *( max(ii[ym]) + min(ii[yp]) ) 
    6775 
    68 print wv 
    69 print sv1, sv2 
    70 print bias 
    71 print bias2 
     76print 'weight vector: ', wv 
     77print 'support vectors: ', sv1, sv2 
     78print 'bias (from points): ', bias 
     79print 'bias (with vectors): ', bias2 
    7280 
    73 # Eqn of hyperplane: 
    74 # wv[0] x + wv[1] y + bias = 0 
    75 hx = array([0, 1.2]) 
     81# plot data 
     82pylab.plot(c1[:,0], c1[:,1], 'bo') 
     83pylab.plot(c2[:,0], c2[:,1], 'yo') 
     84 
     85# plot hyperplane: wv[0] x + wv[1] y + bias = 0 
     86xmin,xmax,ymin,ymax = pylab.axis() 
     87hx = array([floor(xmin-.1), ceil(xmax+.1)]) 
    7688hy = -wv[0]/wv[1] * hx - bias/wv[1] 
     89pylab.plot(hx, hy, 'k-') 
     90#pylab.axis([xmin,xmax,ymin,ymax]) 
     91pylab.show() 
    7792 
    78 pylab.plot(hx, hy, 'k-') 
    79 myshow() 
    80  
    81 # $Id$ 
    82 #  
    8393# end of file 

  • mystic/examples/test_svc2.py

    r826 r827  
    22# 
    33# Author: Patrick Hung (patrickh @caltech) 
     4# Author: Mike McKerns (mmckerns @caltech and @uqfoundation) 
    45# Copyright (c) 1997-2015 California Institute of Technology. 
    56# License: 3-clause BSD.  The full license text is available at: 
    … …  
    78""" 
    89Support Vector Classification. Example 2. 
    9 meristem data. 
     10 
     11using meristem data from data files 
    1012""" 
    1113 
    1214from numpy import * 
    13 from scipy.io import read_array 
    14 import qld, quapro 
    1515import pylab 
    1616from mystic.svctools import * 
    17 import os.path, time 
     17import os.path 
    1818 
    19 def myshow(): 
    20     import Image, tempfile 
    21     name = tempfile.mktemp() 
    22     pylab.savefig(name,dpi=150) 
    23     im = Image.open('%s.png' % name) 
    24     im.show() 
     19# a common objective function for solving a QP problem 
     20# (see http://www.mathworks.com/help/optim/ug/quadprog.html) 
     21def objective(x, H, f): 
     22    return 0.5 * dot(dot(x,H),x) + dot(f,x) 
    2523 
    26 c1 = read_array(os.path.join('DATA','g1x.pts')) 
    27 c2 = read_array(os.path.join('DATA','g2.pts')) 
    28 c1[:,0] += 5 # to make the two sets overlap a little 
     24# SETTINGS 
     25reduced = True  # use a subset of the full data 
     26overlap = False # reduce the distance between the datasets 
    2927 
    30 # interchange ? 
    31 #c1, c2 = c2, c1 
     28c1 = loadtxt(os.path.join('DATA','g1.pts')) 
     29c2 = loadtxt(os.path.join('DATA','g2.pts')) 
     30 
     31if reduced: 
     32    c1 = c1[c1[:,0] > 245] 
     33    c2 = c2[c2[:,0] < 280] 
     34if overlap: c1[:,0] += 3 # make the datasets overlap a little 
    3235 
    3336# the Kernel Matrix (with the linear kernel) 
    34 # Q = multiply.outer(X,X) <--- unfortunately only works when X is a list of scalars... 
    35 # In Mathematica, this would be implemented simply via Outer[K, X, X, 1] 
     37# Q = multiply.outer(X,X)   # NOTE: requires X is a list of scalars 
    3638 
    3739XX = concatenate([c1,-c2]) 
    … …  
    4547f = b 
    4648Aeq = concatenate([ones(c1.shape[0]), -ones(c2.shape[0])]).reshape(1,nx) 
    47 Beq = array([0]) 
     49Beq = array([0.]) 
    4850lb = zeros(nx) 
    49 ub = zeros(nx) + 100000 
     51ub = 1 * ones(nx) 
     52_b = .1 * ones(nx) # good starting value if most solved xi should be zero 
    5053 
    51 #alpha = qld.quadprog2(H, f, None, None, Aeq, Beq, lb, ub) 
    52 alpha,xxx = quapro.quadprog(H, f, None, None, Aeq, Beq, lb, ub) 
    53 #t1 = time.time() 
    54 #for i in range(1000): 
    55 #   alpha,xxx = quapro.quadprog(H, f, None, None, Aeq, Beq, lb, ub) 
    56 #   #alpha = qld.quadprog2(H, f, None, None, Aeq, Beq, lb, ub) 
    57 #t2 = time.time() 
    58 #print "1000 calls to QP took %0.3f s" % (t2-t1) 
    59 print alpha 
     54from mystic.symbolic import linear_symbolic, solve, \ 
     55     generate_solvers as solvers, generate_constraint as constraint 
     56constrain = linear_symbolic(Aeq,Beq) 
     57#NOTE: assumes a single equation of the form: '1.0*x0 + ... = 0.0\n' 
     58x0,rhs = constrain.strip().split(' = ') 
     59x0,xN = x0.split(' + ', 1)  
     60N,x0 = x0.split("*") 
     61constrain = "{x0} = ({rhs} - ({xN}))/{N}".format(x0=x0, xN=xN, N=N, rhs=rhs) 
     62#NOTE: end slight hack (as mystic.symbolic.solve takes __forever__) 
     63constrain = constraint(solvers(constrain)) 
     64 
     65from mystic import supressed 
     66@supressed(5e-2) 
     67def conserve(x): 
     68    return constrain(x) 
     69 
     70from mystic.monitors import VerboseMonitor 
     71mon = VerboseMonitor(10) 
     72 
     73from mystic.solvers import diffev 
     74alpha = diffev(objective, zip(lb,_b), args=(H,f), npop=nx*3, gtol=200,\ 
     75               itermon=mon, \ 
     76               ftol=1e-8, bounds=zip(lb,ub), constraints=conserve, disp=1) 
     77 
     78print 'solved x: ', alpha 
     79print "constraint A*x == 0: ", inner(Aeq, alpha) 
     80print "minimum 0.5*x'Hx + f'*x: ", objective(alpha, H, f) 
    6081 
    6182# the labels and the points 
    … …  
    6485 
    6586wv = WeightVector(alpha, X, y) 
    66 sv1, sv2 = SupportVectors(alpha,y, eps=1e-6) 
     87sv1, sv2 = SupportVectors(alpha,y,eps=1e-6) 
    6788bias = Bias(alpha, X, y) 
    6889 
    69 print wv 
    70 print sv1, sv2 
    71 print bias 
     90ym = (y.flatten()<0).nonzero()[0] 
     91yp = (y.flatten()>0).nonzero()[0] 
     92ii = inner(wv, X) 
     93bias2 = -0.5 *( max(ii[ym]) + min(ii[yp]) ) 
    7294 
    73 # Eqn of hyperplane: 
    74 # wv[0] x + wv[1] y + bias = 0 
     95print 'weight vector: ', wv 
     96print 'support vectors: ', sv1, sv2 
     97print 'bias (from points): ', bias 
     98print 'bias (with vectors): ', bias2 
    7599 
    76 pylab.plot(c1[:,0], c1[:,1], 'ko',markersize=2) 
    77 pylab.plot(c2[:,0], c2[:,1], 'ro',markersize=2) 
     100# plot data 
     101pylab.plot(c1[:,0], c1[:,1], 'bo', markersize=5) 
     102pylab.plot(c2[:,0], c2[:,1], 'yo', markersize=5) 
     103 
     104# plot hyperplane: wv[0] x + wv[1] y + bias = 0 
    78105xmin,xmax,ymin,ymax = pylab.axis() 
    79 hx = array([xmin, xmax]) 
     106hx = array([floor(xmin-.1), ceil(xmax+.1)]) 
    80107hy = -wv[0]/wv[1] * hx - bias/wv[1] 
    81108pylab.plot(hx, hy, 'k-') 
    82 pylab.axis([xmin,xmax,ymin,ymax]) 
    83 pylab.plot(XX[sv1,0], XX[sv1,1], 'ko',markersize=5) 
    84 pylab.plot(-XX[sv2,0], -XX[sv2,1], 'ro',markersize=5) 
    85 pylab.axis('equal') 
    86 myshow() 
     109#pylab.axis([xmin,xmax,ymin,ymax]) 
    87110 
    88 # $Id$ 
    89 #  
     111# plot the support points 
     112pylab.plot(XX[sv1,0], XX[sv1,1], 'b^',markersize=10) 
     113pylab.plot(-XX[sv2,0], -XX[sv2,1], 'y^',markersize=10) 
     114#pylab.axis('equal') 
     115pylab.show() 
     116 
    90117# end of file 
Note: See TracChangeset for help on using the changeset viewer.