+ resolved an issue with the memory mapping in filterReads
[qpalma.git] / qpalma / SIQP.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 import pdb
10
11 import cvxopt.base as cb
12
13 import logging
14 logging.basicConfig(level=logging.DEBUG,format='%(levelname)s %(message)s')
15
16 import qpalma.Configuration as Configuration
17
18 class SIQP:
19 """
20 This class is a wrapper for the cvxopt/cplex qp solvers.
21 It has the purpose to support an iterative addition of
22 constraints to the original problem.
23
24 We want to solve a problem of the form
25
26 min 1/2 * C * x' * P * x + q' * x
27 x
28
29 s.t.
30 < w,f+ > - < w,f- > >= 1 - xi_i for all i, for all f-
31
32 where x is a vector storing information on the parameters w_i and the slacks xi_i
33 so x = [w;xi]
34
35 """
36
37 def __init__(self,fSize,numExamples,c):
38 assert fSize > 0, 'You have to have at least one feature!'
39 assert numExamples > 0, 'You have to have at least one example!'
40 self.numFeatures = fSize
41 self.numExamples = numExamples
42 self.C = c
43 self.EPSILON = 10e-15
44
45 logging.debug("Starting SIQP solver with %d examples. A feature space dim. of %d.", numExamples,fSize)
46 logging.debug("Regularization parameters are C=%f."%(self.C))
47
48 self.old_w = None
49
50 self.dimP = self.numFeatures + self.numExamples
51 self.createRegularizer()
52
53 def createUnitRegularizer(self):
54 # set regularizer to zero
55 self.P = cb.matrix(0.0,(self.dimP,self.dimP))
56 for i in range(self.numFeatures):
57 self.P[i,i] = 1.0
58 # end of zeroing regularizer
59
60 def createSmoothnessRegularizer(self):
61 # set regularizer to zero
62 self.P = cb.matrix(0.0,(self.dimP,self.dimP))
63 for i in range(self.numFeatures):
64 self.P[i,i] = 1.0
65
66 lengthSP = Configuration.numLengthSuppPoints
67 donSP = Configuration.numDonSuppPoints
68 accSP = Configuration.numAccSuppPoints
69 mmatrixSP = Configuration.sizeMatchmatrix[0]\
70 *Configuration.sizeMatchmatrix[1]
71 numq = Configuration.numQualSuppPoints
72 totalQualSP = Configuration.totalQualSuppPoints
73 numQPlifs = Configuration.numQualPlifs
74
75 #lengthGroupParam = 5*1e-9
76 #spliceGroupParam = 1e-9
77
78 lengthGroupParam = 1.0
79 spliceGroupParam = 1.0
80
81 self.P[0:lengthSP,0:lengthSP] *= lengthGroupParam
82
83 beg = lengthSP
84 end = lengthSP+donSP+accSP
85 self.P[beg:end,beg:end] *= spliceGroupParam
86
87 regC = self.numFeatures / (1.0*self.numExamples)
88
89 cfactor = 1.0
90
91 #cfactor = 1e-8
92 #cfactor = 5e-9
93
94 for j in range(lengthSP-1):
95 self.P[j,j+1] = -cfactor
96 self.P[j+1,j] = -cfactor
97 self.P[j,j] += cfactor
98 self.P[j+1,j+1] += cfactor
99
100 #cfactor = 1e-9
101
102 for j in range(lengthSP,lengthSP+donSP-1):
103 cfactor = spliceGroupParam
104 self.P[j,j+1] = -cfactor
105 self.P[j+1,j] = -cfactor
106 self.P[j,j] += cfactor
107 self.P[j+1,j+1] += cfactor
108
109 for j in range(lengthSP+donSP,lengthSP+donSP+accSP-1):
110 cfactor = spliceGroupParam
111 self.P[j,j+1] = -cfactor
112 self.P[j+1,j] = -cfactor
113 self.P[j,j] += cfactor
114 self.P[j+1,j+1] += cfactor
115
116 offset = lengthSP+donSP+accSP+mmatrixSP
117 for k in range(numQPlifs):
118 begin = offset+(k*numq)
119 end = offset+((k+1)*numq)
120
121 for j in range(begin,end-1):
122 self.P[j,j+1] = -1.0
123 self.P[j+1,j] = -1.0
124 self.P[j,j] += 1.0
125
126 # 0.25 for each was already good
127
128 matchGroupParam = 1.0
129 qualityGroupParam = 1.0
130
131 # lengthGroupParam = 0.005
132 # spliceGroupParam = 0.005
133 # matchGroupParam = 0.495
134 # qualityGroupParam = 0.495
135
136 #lengthGroupParam = lengthSP / (1.0*lengthSP+donSP+accSP)
137 #spliceGroupParam = donSP+accSP / (1.0*lengthSP+donSP+accSP)
138 #matchGroupParam = mmatrixSP/(1.0*mmatrixSP+totalQualSP)
139 #qualityGroupParam = totalQualSP/(1.0*mmatrixSP+totalQualSP)
140
141
142 #denom = (1.0*lengthSP+donSP+accSP+mmatrixSP+totalQualSP)
143 #ratio1 = lengthSP+donSP+accSP / denom
144 #ratio2 = mmatrixSP+totalQualSP / denom
145
146 #lengthGroupParam *= ratio1
147 #spliceGroupParam *= ratio1
148
149 #matchGroupParam *= ratio2
150 #qualityGroupParam *= ratio2
151
152 beg = lengthSP+donSP+accSP
153 end = lengthSP+donSP+accSP+mmatrixSP
154 self.P[beg:end,beg:end] *= matchGroupParam
155
156 beg = lengthSP+donSP+accSP+mmatrixSP
157 self.P[beg:-self.numExamples,beg:-self.numExamples] *= qualityGroupParam
158
159 self.P[0:-self.numExamples,0:-self.numExamples] *= 1.0
160
161
162 def createRegularizer(self):
163 #self.createUnitRegularizer()
164 self.createSmoothnessRegularizer()
165
166 q_array = [0.0]*self.numFeatures + [1.0]*self.numExamples
167 self.q = cb.matrix(q_array,(self.numFeatures+self.numExamples,1),'d')
168 self.q = self.C * self.q
169
170 self.G = cb.matrix(0.0,(self.numExamples,self.numFeatures+self.numExamples))
171 for i in range(self.numExamples):
172 self.G[i,self.numFeatures+i] = -1
173
174 self.h = cb.matrix(0,(self.numExamples,1),'d')
175
176 if __name__ == '__main__':
177 sq = SIQP(3,2,100.0)
178
179 """
180 function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
181 % function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
182
183 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184 % qp_solve: min ( <res, f'> + 1/2 res' Q res)
185 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
186 C_qp_penalties = C ; %was 0.01
187 C_qp_smoothness = C ; %was 0.01
188 C_lp_smoothness = 0.00 ; %linear problem parameters
189 C_qp_smallness = 1e-6 ; %was 1e-6
190 column_eps = 1e-3 ;
191
192 Q=zeros(2*126+N) ; % 1/2 * res' * Q * res
193 Q(1:90,1:90) = C_qp_smallness*diag(ones(1,90)) ;
194 Q(91:126,91:126) = C_qp_penalties*diag(ones(1,36)) ;
195 Q(127:2*126,127:2*126) = C_qp_smoothness*diag(ones(1,126)) ;
196 f = [zeros(1,126) C_lp_smoothness*ones(1,126) ones(1,N)/N] ; % <res, f'>
197
198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
199 % qp_solve: LB <= res <= UB
200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
201 INF = 1e20 ;
202 LB = [-INF*ones(1,126) zeros(1,126) zeros(1,N)] ; % lower bound for res
203 UB = [ INF*ones(1,2*126) INF*ones(1,N)] ; % upper bound for res
204
205 """