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