+ resolved some index bugs concerning the intron boundaries resp. splicesite
[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
79 # create CPLEX problem, add objective and constraints to it
80 self.lp = CPX.createprob(self.env, 'test1')
81
82 self.obj = FloatMatrix([0.0]*self.numFeatures + [1.0]*self.numExamples).data
83 cFactor = 1.0*self.C/self.numExamples
84 self.obj = cFactor * self.obj
85 assert self.obj.shape == (self.numVariables,1)
86
87 self.ub = FloatMatrix([self.Inf] * self.numVariables).data
88 self.lb = FloatMatrix([-self.Inf] * self.numFeatures + [0.0]*self.numExamples).data
89
90 self.matbeg = IntMatrix([1]*self.numVariables).data
91 self.matbeg[0] = 0
92
93 self.matcnt = IntMatrix([0]*self.numVariables).data
94 self.matcnt[0] = 1
95
96 self.matind = IntMatrix([0]).data
97 self.matval = FloatMatrix([0.0]).data
98
99 self.rhs = FloatMatrix([0.0]*self.numVariables).data
100 self.sense = CharMatrix(['G']*self.numVariables).data
101
102 self.numConstraints = 1
103
104 CPX.copylp(self.env, self.lp, self.numVariables, self.numConstraints, self.objsen,
105 self.obj, self.rhs, self.sense,
106 self.matbeg, self.matcnt, self.matind, self.matval,
107 self.lb, self.ub)
108
109 self.qmatbeg, self.qmatcnt, self.qmatind, self.qmatval = self.cpx_matrix(self.P)
110 CPX.copyquad(self.env, self.lp, self.qmatbeg, self.qmatcnt, self.qmatind, self.qmatval)
111
112 def addConstraint(self, energy_deltas, example_idx):
113 #def addConstraint(self, energy_deltas, loss, example_idx, useMarginRescaling):
114 """
115 This method adds one contraint to the optimization problem.
116 """
117 self.lastDelta = energy_deltas
118 self.lastIndex = example_idx
119
120 assert -1 < example_idx and example_idx < self.numExamples
121
122 if self.old_w != None:
123 print >> self.protocol, 'Scalar prod. contrib is %f'%((self.old_w[0:self.numFeatures].T*energy_deltas)[0])
124
125 numNewRows = 1
126 numNewCols = 0
127
128 row_rhs = FloatMatrix([1.0]).data
129 row_sense = CharMatrix(['G']).data
130
131 row_matbeg = IntMatrix([0]).data
132
133 row_entries = [0.0]*self.numVariables
134 indices = [0]*self.numVariables
135
136 currentVal = 0.0
137 numNonZeroElems = 0
138 for idx in range(self.numFeatures):
139 currentVal = energy_deltas[idx]
140 if math.fabs(currentVal) > self.EPSILON:
141 row_entries[numNonZeroElems] = currentVal
142 indices[numNonZeroElems] = idx
143 numNonZeroElems += 1
144
145 row_matind = IntMatrix([0]*(numNonZeroElems+1)).data
146 row_matval = FloatMatrix([0.0]*(numNonZeroElems+1)).data
147
148 for pos in range(numNonZeroElems):
149 row_matind[pos] = indices[pos]
150 row_matval[pos] = row_entries[pos]
151
152 # coefficient 1.0 of xi
153 row_matind[numNonZeroElems] = self.numFeatures+example_idx
154 row_matval[numNonZeroElems] = 1.0
155
156 status = CPX.addrows (self.env, self.lp, numNewCols, numNewRows, numNonZeroElems+1, row_rhs,
157 row_sense, row_matbeg, row_matind, row_matval)
158
159 assert status == 0, 'Error ocurred while adding constraint.'
160
161 #CPX.writeprob(self.env,self.lp,"problem_%d.SAV"%self.addedConstraints)
162 #CPX.writeprob(self.env,self.lp,"problem_%d.LP"%self.addedConstraints)
163 self.addedConstraints += 1
164
165 return True
166
167 def solve(self):
168 # solve the problem
169 CPX.qpopt(self.env,self.lp)
170
171 # get the solution
172 lpstat = CPX.getstat(self.env,self.lp)
173 if lpstat in OPTIMAL or lpstat == 6:
174 x, objval = CPX.getx(self.env,self.lp), CPX.getobjval(self.env,self.lp)
175 print >> self.protocol, 'SIQP_CPX: Objective is: %f'%objval
176
177 # Check that objective value is monotonically increasing
178 if math.fabs(objval - self.old_objective_value) >= 10e-5:
179 assert objval >= self.old_objective_value
180
181 self.old_objective_value = objval
182
183 new_w = cb.matrix([0.0]*self.numVariables,(self.numVariables,1))
184 for i in range(self.numVariables):
185 new_w[i,0] = x[i]
186
187 self.old_w = new_w
188
189 # Check for non-negative slacks
190 sumOfSlacks = 0.0
191 print >> self.protocol, 'Slacks are:'
192 for i in range(self.numExamples):
193 print >> self.protocol, '%f ' % new_w[self.numFeatures+i,0]
194 assert new_w[self.numFeatures+i,0] >= 0.0
195 sumOfSlacks += new_w[self.numFeatures+i,0]
196
197 scalarProdDiff = (new_w[0:self.numFeatures,0].T * self.lastDelta)[0]
198 # print >> self.protocol, 'Constraint: %f >= %d - %f'%(scalarProdDiff,self.lastLoss,new_w[self.numFeatures+self.lastIndex])
199
200 # evalute parts of the objective
201 regularizationPart = (0.5 *new_w.T*self.P*new_w)[0]
202 lossPart = 1.0*self.C/self.numExamples * sumOfSlacks
203 assert lossPart >= 0.0, 'Error we have a negative loss contribution'
204 print >> self.protocol, 'Parts of objective: %e %e'%(regularizationPart,lossPart)
205
206 return objval,new_w[0:self.numFeatures,0],new_w[self.numFeatures:,0]
207
208 elif lpstat in UNBOUNDED:
209 print >> self.protocol, "Solution unbounded"
210 assert False
211 elif lpstat in INFEASIBLE:
212 print >> self.protocol, "Solution infeasible"
213 assert False
214 else:
215 print >> self.protocol, "Solution code: ", lpstat
216 assert False
217
218 def cpx_matrix(self,M):
219 """
220 This method creates a sparse CPLEX representation of a given
221 cvxopt Matrix M.
222 """
223
224 nRows,nCols = M.size
225 numEntries = nRows*nCols
226
227 matbeg = IntMatrix([0]*nCols).data
228 matcnt = IntMatrix([0]*nCols).data
229
230 nonZeros = 0
231 elemCounter = 0
232 entry = 0.0
233
234 valVec = [0.0]*numEntries
235 indVec = [0.0]*numEntries
236
237 for currentCol in range(nCols):
238 nonZeros = 0
239 matbeg[currentCol] = elemCounter
240
241 for currentRow in range(nRows):
242 entry = M[currentRow,currentCol]
243 if math.fabs(entry) > self.EPSILON:
244 nonZeros += 1
245 indVec[elemCounter] = currentRow
246 valVec[elemCounter] = entry
247 elemCounter += 1
248
249 matcnt[currentCol] = nonZeros
250
251 totalNonZeros = elemCounter
252
253 matind = IntMatrix([0]*totalNonZeros).data
254 matval = FloatMatrix([0.0]*totalNonZeros).data
255
256 for idx in range(totalNonZeros):
257 matind[idx] = indVec[idx];
258 matval[idx] = valVec[idx];
259
260 print >> self.protocol, 'nRows,nCols: %d,%d'%(nRows,nCols)
261 print >> self.protocol, 'totalNonZeros ' + str(totalNonZeros)
262
263 return matbeg,matcnt,matind,matval
264
265 def cpx_matrix_id(self):
266 """
267 This method creates a sparse CPLEX representation of a given
268 cvxopt Matrix M.
269 """
270
271 nRows,nCols = self.numVariables,self.numVariables
272 matbeg = IntMatrix([0]*nCols).data
273 matcnt = IntMatrix([0]*nCols).data
274
275 matind = IntMatrix([0]*nCols).data
276 matval = FloatMatrix([0.0]*nCols).data
277
278 for i in range(nCols):
279 matbeg[i] = i
280 matcnt[i] = 1
281 matind[i] = i
282 matval[i] = 1.0
283
284 return matbeg,matcnt,matind,matval
285
286 def delete(self):
287 # free the problem
288 CPX.freeprob(self.env,self.lp)
289 # close CPLEX environment
290 CPX.closeCPLEX(self.env)
291
292 def test():
293 fh = open('siqp_cpx.log','w+')
294
295 numFeatures = 3
296 numExamples = 2
297 s = SIQPSolver(numFeatures,numExamples,10.0,fh)
298
299 d1 = cb.matrix([1,0,2],(numFeatures,1))
300 d2 = cb.matrix([1,1,0],(numFeatures,1))
301
302 s.addConstraint(d1,200.0,0,True)
303 s.addConstraint(d2,105.0,1,True)
304
305 w,slacks = s.solve()
306 print w
307 print slacks
308
309 print w.T * d1 + slacks[0]
310 print w.T * d2 + slacks[1]
311
312
313 s.delete()
314 del s
315 fh.close()
316
317 if __name__ == '__main__':
318 test()