+ refactored code further
[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,run):
60 SIQP.__init__(self,fSize,numExamples,c,run)
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 # set intron length to zero
92 #import pdb
93 #pdb.set_trace()
94
95 #self.ub[0] = 0.0
96 #self.ub[1] = 0.0
97
98 #self.lb[0] = 0.0
99 #self.lb[1] = 0.0
100
101 #pdb.set_trace()
102
103 ############################
104
105 self.matbeg = IntMatrix([1]*self.numVariables).data
106 self.matbeg[0] = 0
107
108 self.matcnt = IntMatrix([0]*self.numVariables).data
109 self.matcnt[0] = 1
110
111 self.matind = IntMatrix([0]).data
112 self.matval = FloatMatrix([0.0]).data
113
114 self.rhs = FloatMatrix([0.0]*self.numVariables).data
115 self.sense = CharMatrix(['G']*self.numVariables).data
116
117 self.numConstraints = 1
118
119 CPX.copylp(self.env, self.lp, self.numVariables, self.numConstraints, self.objsen,
120 self.obj, self.rhs, self.sense,
121 self.matbeg, self.matcnt, self.matind, self.matval,
122 self.lb, self.ub)
123
124 import pdb
125 #pdb.set_trace()
126
127 assert sum(self.P[self.numFeatures:,self.numFeatures:]) == 0.0, 'slack variables are regularized'
128
129 self.qmatbeg, self.qmatcnt, self.qmatind, self.qmatval = self.cpx_matrix(self.P)
130 CPX.copyquad(self.env, self.lp, self.qmatbeg, self.qmatcnt, self.qmatind, self.qmatval)
131
132 def enforceMonotonicity(self, begin, end):
133 """
134 This method enforces monotonicity on the parameters indexed by begin to
135 end.
136 """
137 assert -1 < begin < end < self.numFeatures
138
139 for idx in xrange(begin,end):
140 self.enforcePairwiseMonotonicity(idx,idx+1)
141
142 def enforcePairwiseMonotonicity(self,pos1,pos2):
143 energy_deltas = IntMatrix([0.0]*self.numVariables).data
144
145 energy_deltas[pos1] = -1.0
146 energy_deltas[pos2] = 1.0
147
148 numNewRows = 1
149 numNewCols = 0
150
151 row_rhs = FloatMatrix([0.0]).data
152 row_sense = CharMatrix(['G']).data
153
154 row_matbeg = IntMatrix([0]).data
155
156 row_entries = [0.0]*self.numVariables
157 indices = [0]*self.numVariables
158
159 currentVal = 0.0
160 numNonZeroElems = 0
161 for idx in range(self.numFeatures):
162 currentVal = energy_deltas[idx]
163 if math.fabs(currentVal) > self.EPSILON:
164 row_entries[numNonZeroElems] = currentVal
165 indices[numNonZeroElems] = idx
166 numNonZeroElems += 1
167
168 row_matind = IntMatrix([0]*(numNonZeroElems+1)).data
169 row_matval = FloatMatrix([0.0]*(numNonZeroElems+1)).data
170
171 for pos in range(numNonZeroElems):
172 row_matind[pos] = indices[pos]
173 row_matval[pos] = row_entries[pos]
174
175 status = CPX.addrows (self.env, self.lp, numNewCols, numNewRows, numNonZeroElems+1, row_rhs,
176 row_sense, row_matbeg, row_matind, row_matval)
177
178 assert status == 0, 'Error ocurred while adding constraint.'
179
180 self.addedConstraints += 1
181
182 return True
183
184 def addConstraint(self, energy_deltas, example_idx):
185 #def addConstraint(self, energy_deltas, loss, example_idx, useMarginRescaling):
186 """
187 This method adds one contraint to the optimization problem.
188 """
189 self.lastDelta = energy_deltas
190 self.lastIndex = example_idx
191
192 assert -1 < example_idx and example_idx < self.numExamples
193
194 if self.old_w != None:
195 print >> self.protocol, 'Scalar prod. contrib is %f'%((self.old_w[0:self.numFeatures].T*energy_deltas)[0])
196
197 numNewRows = 1
198 numNewCols = 0
199
200 row_rhs = FloatMatrix([1.0]).data
201 row_sense = CharMatrix(['G']).data
202
203 row_matbeg = IntMatrix([0]).data
204
205 row_entries = [0.0]*self.numVariables
206 indices = [0]*self.numVariables
207
208 currentVal = 0.0
209 numNonZeroElems = 0
210 for idx in range(self.numFeatures):
211 currentVal = energy_deltas[idx]
212 if math.fabs(currentVal) > self.EPSILON:
213 row_entries[numNonZeroElems] = currentVal
214 indices[numNonZeroElems] = idx
215 numNonZeroElems += 1
216
217 row_matind = IntMatrix([0]*(numNonZeroElems+1)).data
218 row_matval = FloatMatrix([0.0]*(numNonZeroElems+1)).data
219
220 for pos in range(numNonZeroElems):
221 row_matind[pos] = indices[pos]
222 row_matval[pos] = row_entries[pos]
223
224 # coefficient 1.0 of xi
225 row_matind[numNonZeroElems] = self.numFeatures+example_idx
226 row_matval[numNonZeroElems] = 1.0
227
228 status = CPX.addrows (self.env, self.lp, numNewCols, numNewRows, numNonZeroElems+1, row_rhs,
229 row_sense, row_matbeg, row_matind, row_matval)
230
231 assert status == 0, 'Error ocurred while adding constraint.'
232
233 #CPX.writeprob(self.env,self.lp,"problem_%d.SAV"%self.addedConstraints)
234 #CPX.writeprob(self.env,self.lp,"problem_%d.LP"%self.addedConstraints)
235 self.addedConstraints += 1
236
237 return True
238
239 def solve(self):
240 # solve the problem
241 CPX.qpopt(self.env,self.lp)
242
243 # get the solution
244 lpstat = CPX.getstat(self.env,self.lp)
245 if lpstat in OPTIMAL or lpstat == 6:
246 x, objval = CPX.getx(self.env,self.lp), CPX.getobjval(self.env,self.lp)
247 print >> self.protocol, 'SIQP_CPX: Objective is: %f'%objval
248
249 # Check that objective value is monotonically increasing
250 #if math.fabs(objval - self.old_objective_value) >= 1e-4:
251 # assert objval >= self.old_objective_value
252
253 self.old_objective_value = objval
254
255 new_w = cb.matrix([0.0]*self.numVariables,(self.numVariables,1))
256 for i in range(self.numVariables):
257 new_w[i,0] = x[i]
258
259 self.old_w = new_w
260
261 # Check for non-negative slacks
262 sumOfSlacks = 0.0
263 print >> self.protocol, 'Slacks are:'
264 for i in range(self.numExamples):
265 print >> self.protocol, '%f ' % new_w[self.numFeatures+i,0]
266 assert new_w[self.numFeatures+i,0] >= 0.0
267 sumOfSlacks += new_w[self.numFeatures+i,0]
268
269 scalarProdDiff = (new_w[0:self.numFeatures,0].T * self.lastDelta)[0]
270 # print >> self.protocol, 'Constraint: %f >= %d - %f'%(scalarProdDiff,self.lastLoss,new_w[self.numFeatures+self.lastIndex])
271
272 # evalute parts of the objective
273 regularizationPart = (0.5 *new_w.T*self.P*new_w)[0]
274 lossPart = 1.0*self.C/self.numExamples * sumOfSlacks
275 assert lossPart >= 0.0, 'Error we have a negative loss contribution'
276 print >> self.protocol, 'Parts of objective: %e %e'%(regularizationPart,lossPart)
277
278 return objval,new_w[0:self.numFeatures,0],new_w[self.numFeatures:,0]
279
280 elif lpstat in UNBOUNDED:
281 print >> self.protocol, "Solution unbounded"
282 assert False
283 elif lpstat in INFEASIBLE:
284 print >> self.protocol, "Solution infeasible"
285 assert False
286 else:
287 print >> self.protocol, "Solution code: ", lpstat
288 assert False
289
290 def cpx_matrix(self,M):
291 """
292 This method creates a sparse CPLEX representation of a given
293 cvxopt Matrix M.
294 """
295
296 nRows,nCols = M.size
297 numEntries = nRows*nCols
298
299 matbeg = IntMatrix([0]*nCols).data
300 matcnt = IntMatrix([0]*nCols).data
301
302 nonZeros = 0
303 elemCounter = 0
304 entry = 0.0
305
306 valVec = [0.0]*numEntries
307 indVec = [0.0]*numEntries
308
309 for currentCol in range(nCols):
310 nonZeros = 0
311 matbeg[currentCol] = elemCounter
312
313 for currentRow in range(nRows):
314 entry = M[currentRow,currentCol]
315 if math.fabs(entry) > self.EPSILON:
316 nonZeros += 1
317 indVec[elemCounter] = currentRow
318 valVec[elemCounter] = entry
319 elemCounter += 1
320
321 matcnt[currentCol] = nonZeros
322
323 totalNonZeros = elemCounter
324
325 matind = IntMatrix([0]*totalNonZeros).data
326 matval = FloatMatrix([0.0]*totalNonZeros).data
327
328 for idx in range(totalNonZeros):
329 matind[idx] = indVec[idx];
330 matval[idx] = valVec[idx];
331
332 print >> self.protocol, 'nRows,nCols: %d,%d'%(nRows,nCols)
333 print >> self.protocol, 'totalNonZeros ' + str(totalNonZeros)
334
335 return matbeg,matcnt,matind,matval
336
337 def cpx_matrix_id(self):
338 """
339 This method creates a sparse CPLEX representation of a given
340 cvxopt Matrix M.
341 """
342
343 nRows,nCols = self.numVariables,self.numVariables
344 matbeg = IntMatrix([0]*nCols).data
345 matcnt = IntMatrix([0]*nCols).data
346
347 matind = IntMatrix([0]*nCols).data
348 matval = FloatMatrix([0.0]*nCols).data
349
350 for i in range(nCols):
351 matbeg[i] = i
352 matcnt[i] = 1
353 matind[i] = i
354 matval[i] = 1.0
355
356 return matbeg,matcnt,matind,matval
357
358 def delete(self):
359 # free the problem
360 CPX.freeprob(self.env,self.lp)
361 # close CPLEX environment
362 CPX.closeCPLEX(self.env)
363
364 def test():
365 fh = open('siqp_cpx.log','w+')
366
367 numFeatures = 3
368 numExamples = 2
369 s = SIQPSolver(numFeatures,numExamples,10.0,fh)
370
371 d1 = cb.matrix([1,0,2],(numFeatures,1))
372 d2 = cb.matrix([1,1,0],(numFeatures,1))
373
374 s.addConstraint(d1,200.0,0,True)
375 s.addConstraint(d2,105.0,1,True)
376
377 w,slacks = s.solve()
378 print w
379 print slacks
380
381 print w.T * d1 + slacks[0]
382 print w.T * d2 + slacks[1]
383
384
385 s.delete()
386 del s
387 fh.close()
388
389 if __name__ == '__main__':
390 test()