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