+ extended evaluation function to return positionwise deviations
[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 regC = self.numFeatures / 1.0*self.numExamples
76
77 for j in range(lengthSP-1):
78 self.P[j,j+1] = -1.0
79 self.P[j+1,j] = -1.0
80 self.P[j,j] += 1.0
81
82 for j in range(lengthSP,lengthSP+donSP-1):
83 self.P[j,j+1] = -1.0
84 self.P[j+1,j] = -1.0
85 self.P[j,j] += 1.0
86
87 for j in range(lengthSP+donSP,lengthSP+donSP+accSP-1):
88 self.P[j,j+1] = -1.0
89 self.P[j+1,j] = -1.0
90 self.P[j,j] += 1.0
91
92 offset = lengthSP+donSP+accSP+mmatrixSP
93 for k in range(numQPlifs):
94 begin = offset+(k*numq)
95 end = offset+((k+1)*numq)
96
97 for j in range(begin,end-1):
98 self.P[j,j+1] = -1.0
99 self.P[j+1,j] = -1.0
100 self.P[j,j] += 1.0
101
102 # 0.25 for each was already good
103
104 lengthGroupParam = 0.005
105 spliceGroupParam = 0.005
106
107 matchGroupParam = 0.495
108 qualityGroupParam = 0.495
109
110 self.P[0:lengthSP,0:lengthSP] *= lengthGroupParam
111
112 beg = lengthSP
113 end = lengthSP+donSP+accSP
114 self.P[beg:end,beg:end] *= spliceGroupParam
115
116 beg = lengthSP+donSP+accSP
117 end = lengthSP+donSP+accSP+mmatrixSP
118 self.P[beg:end,beg:end] *= matchGroupParam
119
120 beg = lengthSP+donSP+accSP+mmatrixSP
121 self.P[beg:-self.numExamples,beg:-self.numExamples] *= qualityGroupParam
122
123 self.P[0:-self.numExamples,0:-self.numExamples] *= regC*1e-3
124
125
126 def createRegularizer(self):
127 #self.createUnitRegularizer()
128 self.createSmoothnessRegularizer()
129
130 q_array = [0.0]*self.numFeatures + [1.0]*self.numExamples
131 self.q = cb.matrix(q_array,(self.numFeatures+self.numExamples,1),'d')
132 self.q = self.C * self.q
133
134 self.G = cb.matrix(0.0,(self.numExamples,self.numFeatures+self.numExamples))
135 for i in range(self.numExamples):
136 self.G[i,self.numFeatures+i] = -1
137
138 self.h = cb.matrix(0,(self.numExamples,1),'d')
139
140 if __name__ == '__main__':
141 sq = SIQP(3,2,100.0)
142
143 """
144 function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
145 % function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
146
147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
148 % qp_solve: min ( <res, f'> + 1/2 res' Q res)
149 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
150 C_qp_penalties = C ; %was 0.01
151 C_qp_smoothness = C ; %was 0.01
152 C_lp_smoothness = 0.00 ; %linear problem parameters
153 C_qp_smallness = 1e-6 ; %was 1e-6
154 column_eps = 1e-3 ;
155
156 Q=zeros(2*126+N) ; % 1/2 * res' * Q * res
157 Q(1:90,1:90) = C_qp_smallness*diag(ones(1,90)) ;
158 Q(91:126,91:126) = C_qp_penalties*diag(ones(1,36)) ;
159 Q(127:2*126,127:2*126) = C_qp_smoothness*diag(ones(1,126)) ;
160 f = [zeros(1,126) C_lp_smoothness*ones(1,126) ones(1,N)/N] ; % <res, f'>
161
162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
163 % qp_solve: LB <= res <= UB
164 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
165 INF = 1e20 ;
166 LB = [-INF*ones(1,126) zeros(1,126) zeros(1,N)] ; % lower bound for res
167 UB = [ INF*ones(1,2*126) INF*ones(1,N)] ; % upper bound for res
168
169 """