+ changed check for intial setting of old_w. (matrix != None does not work in cvxopt)
[qpalma.git] / qpalma / SIQP_CVXOPT.py
1 # This program is free software; you can redistribute it and/or modify
2 # it under the terms of the GNU General Public License as published by
3 # the Free Software Foundation; either version 2 of the License, or
4 # (at your option) any later version.
5 #
6 # Written (W) 2008 Fabio De Bona
7 # Copyright (C) 2008 Max-Planck-Society
8
9 import sys
10 import os
11 import os.path
12 import logging
13 import math
14 import pdb
15
16 import cvxopt.base as cb
17 import cvxopt.solvers as cs
18
19 from SIQP import SIQP
20
21 acc = 1e-7
22
23 if __debug__:
24 cs.options['show_progress'] = True # default: True
25 else:
26 cs.options['show_progress'] = False # default: True
27
28 cs.options['maxiters'] = 200 # default: 100
29 cs.options['abstol'] = acc # default: 1e-7
30 cs.options['reltol'] = acc # default: 1e-7
31 cs.options['feastol'] = acc # default: 1e-7
32 #cs.options['refinement'] = True # default: True
33
34
35
36 class SIQPSolver(SIQP):
37 """
38 This class is a wrapper for the cvxopt qp solver.
39 It has the purpose to support an iterative addition of
40 constraints to the original problem.
41
42 """
43
44 def __init__(self,fSize,numExamples,c,logfh,settings):
45 SIQP.__init__(self,fSize,numExamples,c,settings)
46
47 self.numConstraints = 0
48 self.solver_runs = 0
49 self.old_objective_value = -(sys.maxint-1)
50 self.logfh = logfh
51
52 self.old_w = cb.matrix([0.0] ,(1,1))
53
54 #self.P.tofile(open('matrix_P_%d_%d'%(self.P.size[0],self.P.size[1]),'w+'))
55 #self.q.tofile(open('matrix_q_%d_%d'%(self.q.size[0],self.q.size[1]),'w+'))
56 #self.G.tofile(open('matrix_G_%d_%d'%(self.G.size[0],self.G.size[1]),'w+'))
57 #self.h.tofile(open('matrix_h_%d_%d'%(self.h.size[0],self.h.size[1]),'w+'))
58
59 def plog(self,string):
60 if self.logfh != None:
61 self.logfh.write(string)
62 self.logfh.flush()
63 else:
64 print string
65
66 def addConstraint(self, energy_deltas, idx):
67 energy_deltas = cb.matrix(energy_deltas)
68 loss = 1.0
69 useMarginRescaling = True
70 self.currentIdx = idx
71 self.currentLoss = loss
72 self.energy_deltas = cb.matrix(energy_deltas)
73 scalar_prod = 0.0
74
75 energy_deltas = energy_deltas.T
76
77 assert useMarginRescaling
78
79 if loss < self.EPSILON:
80 return False
81
82 if self.old_w.size != (1,1):
83 scalar_prod = energy_deltas * self.old_w[0:self.numFeatures,0]
84 old_slack = self.old_w[self.numFeatures+idx,0]
85 if useMarginRescaling:
86 if scalar_prod[0] + old_slack > loss: # leave some room for inaccuracies
87 print 'Current example already fulfills constraint.'
88 self.plog('Current example already fulfills constraint.\n')
89 return False
90 else:
91 self.plog("constraint at current solution: %f <= -1+%f/%f = %f" %(scalar_prod[0], old_slack, loss, -1 + old_slack/loss) )
92 if scalar_prod[0] < -1+old_slack/loss + 1e-6: # leave some room for inaccuracies
93 self.plog('Current example already fulfills constraint.\n')
94 return False
95
96 rows,cols = self.G.size
97 self.newG = cb.matrix(0.0,(rows+1,cols))
98 self.newG[0:rows,:] = self.G[0:rows,:]
99
100 self.newG[-1,0:self.numFeatures] = -1.0 * energy_deltas
101
102 if useMarginRescaling:
103 # margin rescaling
104 self.newG[-1,self.numFeatures+idx] = -1.0
105 else:
106 # slack rescaling
107 self.newG[-1,self.numFeatures+idx] = -1.0 / loss
108 self.G = self.newG
109
110 rows,cols = self.h.size
111 self.new_h = cb.matrix(0.0,(rows+1,cols))
112 self.new_h[0:rows,0] = self.h[0:rows,0]
113
114 if useMarginRescaling:
115 # margin rescaling
116 self.new_h[-1,0] = -1.0 * loss
117 else:
118 # slack rescaling
119 self.new_h[-1,0] = -1.0
120 self.h = self.new_h
121
122 self.numConstraints = self.G.size[0]
123 self.plog('number of constraints is %d\n' % self.numConstraints)
124
125 assert self.G.size[0] == self.h.size[0]
126
127 return True
128
129 def solve(self):
130 d = cs.qp(self.P,self.q,self.G,self.h)
131 self.solver_runs += 1
132
133 self.plog('Return status of solver: ' + d['status'])
134 new_w = d['x']
135 self.old_w = new_w
136 self.plog('slacks are:\n')
137
138 slacks_string = ''
139 for i in range(self.numFeatures,new_w.size[0]):
140 slacks_string += ('%e ' % new_w[i])
141 assert new_w[i] >= 0, 'Error found a negative slack variable'
142
143 self.plog(slacks_string+'\n')
144 self.plog('end of slacks\n')
145
146 obj_value = 0.5 * new_w.T*self.P*new_w + self.q.T*new_w
147 obj_value = obj_value[0]
148
149 # evalute parts of the objective
150 regularizationPart = (0.5 *new_w.T*self.P*new_w)[0]
151 lossPart = (self.q.T*new_w)[0]
152 assert lossPart >= 0.0, 'Error we have a negative loss contribution'
153 self.plog('Parts of objective: %e %e\n'%(regularizationPart,lossPart))
154
155 print 'SIQP_CVXOPT: Objective is: %e'%obj_value
156 self.plog('SIQP_CVXOPT: Objective is: %e\n'%obj_value)
157
158 print self.old_objective_value
159 print obj_value
160 assert self.old_objective_value <= obj_value + self.EPSILON
161
162 self.old_objective_value = obj_value
163
164 #scalar_prod = self.currentLoss + (self.energy_deltas * new_w[0:self.numFeatures,0])
165 #assert ( scalar_prod <= new_w[self.numFeatures+self.currentIdx] )
166
167 return obj_value,new_w[0:self.numFeatures,0],new_w[self.numFeatures:,0]