d3be2e532abc1f5fe242677ce4a407558476398b
[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 import logging
20 logging.basicConfig(level=logging.DEBUG,format='%(levelname)s %(message)s')
21
22 from SIQP import SIQP
23
24 acc = 1e-7
25
26 if __debug__:
27 cs.options['show_progress'] = True # default: True
28 else:
29 cs.options['show_progress'] = False # default: True
30
31 cs.options['maxiters'] = 200 # default: 100
32 cs.options['abstol'] = acc # default: 1e-7
33 cs.options['reltol'] = acc # default: 1e-7
34 cs.options['feastol'] = acc # default: 1e-7
35 #cs.options['refinement'] = True # default: True
36
37 class SIQPSolver(SIQP):
38 """
39 This class is a wrapper for the cvxopt qp solver.
40 It has the purpose to support an iterative addition of
41 constraints to the original problem.
42
43 """
44
45 def __init__(self,fSize,numExamples,c,proto,run):
46 SIQP.__init__(self,fSize,numExamples,c,run)
47 self.numConstraints = 0
48 self.solver_runs = 0
49 self.old_objective_value = -(sys.maxint-1)
50 self.protocol = proto
51
52 #self.P.tofile(open('matrix_P_%d_%d'%(self.P.size[0],self.P.size[1]),'w+'))
53 #self.q.tofile(open('matrix_q_%d_%d'%(self.q.size[0],self.q.size[1]),'w+'))
54 #self.G.tofile(open('matrix_G_%d_%d'%(self.G.size[0],self.G.size[1]),'w+'))
55 #self.h.tofile(open('matrix_h_%d_%d'%(self.h.size[0],self.h.size[1]),'w+'))
56
57 def addConstraint(self, energy_deltas, idx):
58 energy_deltas = cb.matrix(energy_deltas)
59 loss = 1.0
60 useMarginRescaling = True
61 self.currentIdx = idx
62 self.currentLoss = loss
63 self.energy_deltas = cb.matrix(energy_deltas)
64 scalar_prod = 0.0
65
66 energy_deltas = energy_deltas.T
67
68 assert useMarginRescaling
69
70 if loss < self.EPSILON:
71 return False
72
73 if self.old_w != None:
74 scalar_prod = energy_deltas * self.old_w[0:self.numFeatures,0]
75 old_slack = self.old_w[self.numFeatures+idx,0]
76 if useMarginRescaling:
77 if scalar_prod[0] + old_slack > loss: # leave some room for inaccuracies
78 print 'Current example already fulfills constraint.'
79 print >> self.protocol, 'Current example already fulfills constraint.'
80 return False
81 else:
82 print >> self.protocol, "constraint at current solution: %f <= -1+%f/%f = %f" %(scalar_prod[0], old_slack, loss, -1 + old_slack/loss)
83 if scalar_prod[0] < -1+old_slack/loss + 1e-6: # leave some room for inaccuracies
84 print >> self.protocol, 'Current example already fulfills constraint.'
85 logging.debug('Current example already fulfills constraint.')
86 return False
87
88 rows,cols = self.G.size
89 self.newG = cb.matrix(0.0,(rows+1,cols))
90 self.newG[0:rows,:] = self.G[0:rows,:]
91
92 self.newG[-1,0:self.numFeatures] = -1.0 * energy_deltas
93
94 if useMarginRescaling:
95 # margin rescaling
96 self.newG[-1,self.numFeatures+idx] = -1.0
97 else:
98 # slack rescaling
99 self.newG[-1,self.numFeatures+idx] = -1.0 / loss
100 self.G = self.newG
101
102 rows,cols = self.h.size
103 self.new_h = cb.matrix(0.0,(rows+1,cols))
104 self.new_h[0:rows,0] = self.h[0:rows,0]
105
106 if useMarginRescaling:
107 # margin rescaling
108 self.new_h[-1,0] = -1.0 * loss
109 else:
110 # slack rescaling
111 self.new_h[-1,0] = -1.0
112 self.h = self.new_h
113
114 self.numConstraints = self.G.size[0]
115 print >> self.protocol, 'number of constraints is %d' % self.numConstraints
116
117 assert self.G.size[0] == self.h.size[0]
118
119 return True
120
121 def solve(self):
122 d = cs.qp(self.P,self.q,self.G,self.h)
123 self.solver_runs += 1
124
125 print >> self.protocol, 'Return status of solver: ' + d['status']
126 new_w = d['x']
127 self.old_w = new_w
128 print >> self.protocol, 'slacks are'
129 for i in range(self.numFeatures,new_w.size[0]):
130 print >> self.protocol, "%e " % new_w[i],
131 assert new_w[i] >= 0, 'Error found a negative slack variable'
132 print >> self.protocol, "end of slacks"
133
134 obj_value = 0.5 * new_w.T*self.P*new_w + self.q.T*new_w
135 obj_value = obj_value[0]
136
137 # evalute parts of the objective
138 regularizationPart = (0.5 *new_w.T*self.P*new_w)[0]
139 lossPart = (self.q.T*new_w)[0]
140 assert lossPart >= 0.0, 'Error we have a negative loss contribution'
141 print >> self.protocol, 'Parts of objective: %e %e'%(regularizationPart,lossPart)
142
143 print 'SIQP_CVXOPT: Objective is: %e'%obj_value
144 print >> self.protocol, 'SIQP_CVXOPT: Objective is: %e'%obj_value
145
146 print self.old_objective_value
147 print obj_value
148 assert self.old_objective_value <= obj_value + self.EPSILON
149
150 self.old_objective_value = obj_value
151
152 #scalar_prod = self.currentLoss + (self.energy_deltas * new_w[0:self.numFeatures,0])
153 #assert ( scalar_prod <= new_w[self.numFeatures+self.currentIdx] )
154
155 return obj_value,new_w[0:self.numFeatures,0],new_w[self.numFeatures:,0]
156
157 def test():
158 fh = open('siqp_cvxopt.log','w+')
159 numFeatures = 3
160 numExamples = 2
161 s = SIQPSolver(numFeatures,numExamples,100.0,fh)
162
163 d1 = cb.matrix([1,0,2],(numFeatures,1))
164 d2 = cb.matrix([1,1,0],(numFeatures,1))
165
166 s.addConstraint(d1,10.0,0,True)
167 s.addConstraint(d2,105.0,1,True)
168
169 w,slacks = s.solve()
170 print w
171 print slacks
172
173 print w.T * d1 + slacks[0]
174 print w.T * d2 + slacks[1]
175
176 del s
177 fh.close()
178
179 if __name__ == '__main__':
180 test()