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