+ rearranged code
[qpalma.git] / qpalma / SIQP_CPX.py
1 #!/usr/bin/env python
2 # -*- coding:latin-1 -*-
3
4 import sys
5 import os
6 import os.path
7 import logging
8 import math
9
10 import cvxopt.base as cb
11
12 import logging
13 logging.basicConfig(level=logging.DEBUG,format='%(levelname)s %(message)s')
14
15 import CPX
16 import numpy as ny
17 import numpy.matlib as nyml
18 from cplex_const import *
19
20 from SIQP import SIQP
21
22 OPTIMAL = [CPX_STAT_OPTIMAL, CPXMIP_OPTIMAL, CPXMIP_OPTIMAL_TOL]
23 UNBOUNDED = [CPX_STAT_UNBOUNDED, CPXMIP_UNBOUNDED]
24 INFEASIBLE = [CPX_STAT_INFEASIBLE, CPXMIP_INFEASIBLE, CPXMIP_INForUNBD]
25
26 class FloatMatrix:
27
28 TYPE = nyml.float64
29
30 def __init__(self,list):
31 self.data = nyml.mat(list, dtype=self.TYPE).T
32
33 class IntMatrix:
34
35 TYPE = nyml.int32
36
37 def __init__(self,list):
38 self.data = nyml.mat(list, dtype=self.TYPE).T
39
40 class CharMatrix:
41
42 TYPE = '|S1'
43
44 def __init__(self,list):
45 self.data = nyml.mat(list, dtype=self.TYPE).T
46
47 class SIQPSolver(SIQP):
48 """
49 This algorithm solves the problem
50
51 min 1/2 x'Px + q'x
52 s.t.
53 Gx >= h
54 """
55
56 # Define types and constants
57 Inf = CPX_INFBOUND
58
59 def __init__(self,fSize,numExamples,c,proto):
60 SIQP.__init__(self,fSize,numExamples,c)
61 self.numFeatures = fSize
62 self.numExamples = numExamples
63 self.numVariables = self.numFeatures + self.numExamples
64 self.old_w = None
65 self.old_objective_value = -(sys.maxint-1)
66
67 self.addedConstraints = 0
68
69 self.protocol = proto
70
71 self.objsen = CPX_MIN
72
73 # open CPLEX environment, set some parameters
74 self.env = CPX.openCPLEX()
75
76 CPX.setintparam(self.env, CPX_PARAM_SCRIND, CPX_ON) # print >> self.protocol, info to screen
77 CPX.setintparam(self.env, CPX_PARAM_DATACHECK, CPX_ON)
78 #CPX.setintparam(self.env, CPX_PARAM_QPMETHOD, 2)
79
80 # create CPLEX problem, add objective and constraints to it
81 self.lp = CPX.createprob(self.env, 'test1')
82
83 self.obj = FloatMatrix([0.0]*self.numFeatures + [1.0]*self.numExamples).data
84 cFactor = 1.0*self.C/self.numExamples
85 self.obj = cFactor * self.obj
86 assert self.obj.shape == (self.numVariables,1)
87
88 self.ub = FloatMatrix([self.Inf] * self.numVariables).data
89 self.lb = FloatMatrix([-self.Inf] * self.numFeatures + [0.0]*self.numExamples).data
90
91 self.matbeg = IntMatrix([1]*self.numVariables).data
92 self.matbeg[0] = 0
93
94 self.matcnt = IntMatrix([0]*self.numVariables).data
95 self.matcnt[0] = 1
96
97 self.matind = IntMatrix([0]).data
98 self.matval = FloatMatrix([0.0]).data
99
100 self.rhs = FloatMatrix([0.0]*self.numVariables).data
101 self.sense = CharMatrix(['G']*self.numVariables).data
102
103 self.numConstraints = 1
104
105 CPX.copylp(self.env, self.lp, self.numVariables, self.numConstraints, self.objsen,
106 self.obj, self.rhs, self.sense,
107 self.matbeg, self.matcnt, self.matind, self.matval,
108 self.lb, self.ub)
109
110 import pdb
111 #pdb.set_trace()
112
113 assert sum(self.P[self.numFeatures:,self.numFeatures:]) == 0.0, 'slack variables are regularized'
114
115 self.qmatbeg, self.qmatcnt, self.qmatind, self.qmatval = self.cpx_matrix(self.P)
116 CPX.copyquad(self.env, self.lp, self.qmatbeg, self.qmatcnt, self.qmatind, self.qmatval)
117
118 def enforceMonotonicity(self, begin, end):
119 """
120 This method enforces monotonicity on the parameters indexed by begin to
121 end.
122 """
123 assert -1 < begin < end < self.numFeatures
124
125 for idx in xrange(begin,end):
126 self.enforcePairwiseMonotonicity(idx,idx+1)
127
128 def enforcePairwiseMonotonicity(self,pos1,pos2):
129 energy_deltas = IntMatrix([0.0]*self.numVariables).data
130
131 energy_deltas[pos1] = -1.0
132 energy_deltas[pos2] = 1.0
133
134 numNewRows = 1
135 numNewCols = 0
136
137 row_rhs = FloatMatrix([0.0]).data
138 row_sense = CharMatrix(['G']).data
139
140 row_matbeg = IntMatrix([0]).data
141
142 row_entries = [0.0]*self.numVariables
143 indices = [0]*self.numVariables
144
145 currentVal = 0.0
146 numNonZeroElems = 0
147 for idx in range(self.numFeatures):
148 currentVal = energy_deltas[idx]
149 if math.fabs(currentVal) > self.EPSILON:
150 row_entries[numNonZeroElems] = currentVal
151 indices[numNonZeroElems] = idx
152 numNonZeroElems += 1
153
154 row_matind = IntMatrix([0]*(numNonZeroElems+1)).data
155 row_matval = FloatMatrix([0.0]*(numNonZeroElems+1)).data
156
157 for pos in range(numNonZeroElems):
158 row_matind[pos] = indices[pos]
159 row_matval[pos] = row_entries[pos]
160
161 status = CPX.addrows (self.env, self.lp, numNewCols, numNewRows, numNonZeroElems+1, row_rhs,
162 row_sense, row_matbeg, row_matind, row_matval)
163
164 assert status == 0, 'Error ocurred while adding constraint.'
165
166 self.addedConstraints += 1
167
168 return True
169
170 def addConstraint(self, energy_deltas, example_idx):
171 #def addConstraint(self, energy_deltas, loss, example_idx, useMarginRescaling):
172 """
173 This method adds one contraint to the optimization problem.
174 """
175 self.lastDelta = energy_deltas
176 self.lastIndex = example_idx
177
178 assert -1 < example_idx and example_idx < self.numExamples
179
180 if self.old_w != None:
181 print >> self.protocol, 'Scalar prod. contrib is %f'%((self.old_w[0:self.numFeatures].T*energy_deltas)[0])
182
183 numNewRows = 1
184 numNewCols = 0
185
186 row_rhs = FloatMatrix([1.0]).data
187 row_sense = CharMatrix(['G']).data
188
189 row_matbeg = IntMatrix([0]).data
190
191 row_entries = [0.0]*self.numVariables
192 indices = [0]*self.numVariables
193
194 currentVal = 0.0
195 numNonZeroElems = 0
196 for idx in range(self.numFeatures):
197 currentVal = energy_deltas[idx]
198 if math.fabs(currentVal) > self.EPSILON:
199 row_entries[numNonZeroElems] = currentVal
200 indices[numNonZeroElems] = idx
201 numNonZeroElems += 1
202
203 row_matind = IntMatrix([0]*(numNonZeroElems+1)).data
204 row_matval = FloatMatrix([0.0]*(numNonZeroElems+1)).data
205
206 for pos in range(numNonZeroElems):
207 row_matind[pos] = indices[pos]
208 row_matval[pos] = row_entries[pos]
209
210 # coefficient 1.0 of xi
211 row_matind[numNonZeroElems] = self.numFeatures+example_idx
212 row_matval[numNonZeroElems] = 1.0
213
214 status = CPX.addrows (self.env, self.lp, numNewCols, numNewRows, numNonZeroElems+1, row_rhs,
215 row_sense, row_matbeg, row_matind, row_matval)
216
217 assert status == 0, 'Error ocurred while adding constraint.'
218
219 #CPX.writeprob(self.env,self.lp,"problem_%d.SAV"%self.addedConstraints)
220 #CPX.writeprob(self.env,self.lp,"problem_%d.LP"%self.addedConstraints)
221 self.addedConstraints += 1
222
223 return True
224
225 def solve(self):
226 # solve the problem
227 CPX.qpopt(self.env,self.lp)
228
229 # get the solution
230 lpstat = CPX.getstat(self.env,self.lp)
231 if lpstat in OPTIMAL or lpstat == 6:
232 x, objval = CPX.getx(self.env,self.lp), CPX.getobjval(self.env,self.lp)
233 print >> self.protocol, 'SIQP_CPX: Objective is: %f'%objval
234
235 # Check that objective value is monotonically increasing
236 if math.fabs(objval - self.old_objective_value) >= 10e-5:
237 assert objval >= self.old_objective_value
238
239 self.old_objective_value = objval
240
241 new_w = cb.matrix([0.0]*self.numVariables,(self.numVariables,1))
242 for i in range(self.numVariables):
243 new_w[i,0] = x[i]
244
245 self.old_w = new_w
246
247 # Check for non-negative slacks
248 sumOfSlacks = 0.0
249 print >> self.protocol, 'Slacks are:'
250 for i in range(self.numExamples):
251 print >> self.protocol, '%f ' % new_w[self.numFeatures+i,0]
252 assert new_w[self.numFeatures+i,0] >= 0.0
253 sumOfSlacks += new_w[self.numFeatures+i,0]
254
255 scalarProdDiff = (new_w[0:self.numFeatures,0].T * self.lastDelta)[0]
256 # print >> self.protocol, 'Constraint: %f >= %d - %f'%(scalarProdDiff,self.lastLoss,new_w[self.numFeatures+self.lastIndex])
257
258 # evalute parts of the objective
259 regularizationPart = (0.5 *new_w.T*self.P*new_w)[0]
260 lossPart = 1.0*self.C/self.numExamples * sumOfSlacks
261 assert lossPart >= 0.0, 'Error we have a negative loss contribution'
262 print >> self.protocol, 'Parts of objective: %e %e'%(regularizationPart,lossPart)
263
264 return objval,new_w[0:self.numFeatures,0],new_w[self.numFeatures:,0]
265
266 elif lpstat in UNBOUNDED:
267 print >> self.protocol, "Solution unbounded"
268 assert False
269 elif lpstat in INFEASIBLE:
270 print >> self.protocol, "Solution infeasible"
271 assert False
272 else:
273 print >> self.protocol, "Solution code: ", lpstat
274 assert False
275
276 def cpx_matrix(self,M):
277 """
278 This method creates a sparse CPLEX representation of a given
279 cvxopt Matrix M.
280 """
281
282 nRows,nCols = M.size
283 numEntries = nRows*nCols
284
285 matbeg = IntMatrix([0]*nCols).data
286 matcnt = IntMatrix([0]*nCols).data
287
288 nonZeros = 0
289 elemCounter = 0
290 entry = 0.0
291
292 valVec = [0.0]*numEntries
293 indVec = [0.0]*numEntries
294
295 for currentCol in range(nCols):
296 nonZeros = 0
297 matbeg[currentCol] = elemCounter
298
299 for currentRow in range(nRows):
300 entry = M[currentRow,currentCol]
301 if math.fabs(entry) > self.EPSILON:
302 nonZeros += 1
303 indVec[elemCounter] = currentRow
304 valVec[elemCounter] = entry
305 elemCounter += 1
306
307 matcnt[currentCol] = nonZeros
308
309 totalNonZeros = elemCounter
310
311 matind = IntMatrix([0]*totalNonZeros).data
312 matval = FloatMatrix([0.0]*totalNonZeros).data
313
314 for idx in range(totalNonZeros):
315 matind[idx] = indVec[idx];
316 matval[idx] = valVec[idx];
317
318 print >> self.protocol, 'nRows,nCols: %d,%d'%(nRows,nCols)
319 print >> self.protocol, 'totalNonZeros ' + str(totalNonZeros)
320
321 return matbeg,matcnt,matind,matval
322
323 def cpx_matrix_id(self):
324 """
325 This method creates a sparse CPLEX representation of a given
326 cvxopt Matrix M.
327 """
328
329 nRows,nCols = self.numVariables,self.numVariables
330 matbeg = IntMatrix([0]*nCols).data
331 matcnt = IntMatrix([0]*nCols).data
332
333 matind = IntMatrix([0]*nCols).data
334 matval = FloatMatrix([0.0]*nCols).data
335
336 for i in range(nCols):
337 matbeg[i] = i
338 matcnt[i] = 1
339 matind[i] = i
340 matval[i] = 1.0
341
342 return matbeg,matcnt,matind,matval
343
344 def delete(self):
345 # free the problem
346 CPX.freeprob(self.env,self.lp)
347 # close CPLEX environment
348 CPX.closeCPLEX(self.env)
349
350 def test():
351 fh = open('siqp_cpx.log','w+')
352
353 numFeatures = 3
354 numExamples = 2
355 s = SIQPSolver(numFeatures,numExamples,10.0,fh)
356
357 d1 = cb.matrix([1,0,2],(numFeatures,1))
358 d2 = cb.matrix([1,1,0],(numFeatures,1))
359
360 s.addConstraint(d1,200.0,0,True)
361 s.addConstraint(d2,105.0,1,True)
362
363 w,slacks = s.solve()
364 print w
365 print slacks
366
367 print w.T * d1 + slacks[0]
368 print w.T * d2 + slacks[1]
369
370
371 s.delete()
372 del s
373 fh.close()
374
375 if __name__ == '__main__':
376 test()