+ renaming business
[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 def asymCoords(t):
19 size_upstream = t[0]
20 size_downstream = t[1]
21 assert size_upstream >= 1 and size_upstream <= 39
22 assert size_downstream >= 1 and size_downstream <= 39
23
24 if size_downstream > size_upstream:
25 asym_small = size_upstream
26 asym_big = size_downstream
27 else:
28 asym_small = size_downstream
29 asym_big = size_upstream
30
31 asym_offset = 0
32 for r in range(asym_small):
33 asym_offset += (40-r)-r;
34
35 return (asym_offset+asym_big)-1
36
37 def calcTuples(i,j):
38 return filter(lambda t: t[0]<=t[1] and t[0]>0 and t[0]<40 and t[1]>0 and t[1]<40,[(i+1,j),(i-1,j),(i,j+1),(i,j-1)])
39
40 class SIQP:
41 """
42 This class is a wrapper for the cvxopt qp solver.
43 It has the purpose to support an iterative addition of
44 constraints to the original problem.
45
46 We want to solve a problem of the form
47
48 min 1/2 * C * x' * P * x + q' * x
49 x
50
51 s.t.
52 < w,f+ > - < w,f- > >= 1 - xi_i for all i, for all f-
53
54 where x is a vector storing information on the parameters w_i and the slacks xi_i
55 so x = [w;xi]
56
57 """
58
59 def __init__(self,fSize,numExamples,c):
60 assert fSize > 0, 'You have to have at least one feature!'
61 assert numExamples > 0, 'You have to have at least one example!'
62 self.numFeatures = fSize
63 self.numExamples = numExamples
64 self.C = c
65 self.EPSILON = 10e-15
66
67 logging.debug("Starting SIQP solver with %d examples. A feature space dim. of %d.", numExamples,fSize)
68 logging.debug("Regularization parameters are C=%f."%(self.C))
69
70 self.old_w = None
71
72 self.dimP = self.numFeatures + self.numExamples
73 self.createRegularizer()
74
75 def createUnitRegularizer(self):
76 # set regularizer to zero
77 self.P = cb.matrix(0.0,(self.dimP,self.dimP))
78 for i in range(self.numFeatures):
79 self.P[i,i] = 1.0
80 # end of zeroing regularizer
81
82 def createSmoothnessRegularizer(self):
83 # set regularizer to zero
84 self.P = cb.matrix(0.0,(self.dimP,self.dimP))
85 for i in range(self.numFeatures):
86 self.P[i,i] = 1.0
87
88 lengthSP = Configuration.numLengthSuppPoints
89 donSP = Configuration.numDonSuppPoints
90 accSP = Configuration.numAccSuppPoints
91 mmatrixSP = Configuration.sizeMatchmatrix[0]\
92 *Configuration.sizeMatchmatrix[1]
93 numq = Configuration.numQualSuppPoints
94 totalQualSP = Configuration.totalQualSuppPoints
95 numQPlifs = Configuration.numQualPlifs
96
97 for j in range(lengthSP-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 for j in range(lengthSP,lengthSP+donSP-1):
103 self.P[j,j+1] = -1.0
104 self.P[j+1,j] = -1.0
105 self.P[j,j] += 1.0
106
107 for j in range(lengthSP+donSP,lengthSP+donSP+accSP-1):
108 self.P[j,j+1] = -1.0
109 self.P[j+1,j] = -1.0
110 self.P[j,j] += 1.0
111
112 #for k in range(numQPlifs):
113 #offset = lengthSP+donSP+accSP+mmatrixSP
114 #for j in range(offset+,):
115
116 def createRegularizer(self):
117 #self.createUnitRegularizer()
118 self.createSmoothnessRegularizer()
119
120 q_array = [0.0]*self.numFeatures + [1.0]*self.numExamples
121 self.q = cb.matrix(q_array,(self.numFeatures+self.numExamples,1),'d')
122 self.q = self.C * self.q
123
124 self.G = cb.matrix(0.0,(self.numExamples,self.numFeatures+self.numExamples))
125 for i in range(self.numExamples):
126 self.G[i,self.numFeatures+i] = -1
127
128 self.h = cb.matrix(0,(self.numExamples,1),'d')
129
130 if __name__ == '__main__':
131 sq = SIQP(3,2,100.0)
132
133 """
134 function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
135 % function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
136
137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
138 % qp_solve: min ( <res, f'> + 1/2 res' Q res)
139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140 C_qp_penalties = C ; %was 0.01
141 C_qp_smoothness = C ; %was 0.01
142 C_lp_smoothness = 0.00 ; %linear problem parameters
143 C_qp_smallness = 1e-6 ; %was 1e-6
144 column_eps = 1e-3 ;
145
146 Q=zeros(2*126+N) ; % 1/2 * res' * Q * res
147 Q(1:90,1:90) = C_qp_smallness*diag(ones(1,90)) ;
148 Q(91:126,91:126) = C_qp_penalties*diag(ones(1,36)) ;
149 Q(127:2*126,127:2*126) = C_qp_smoothness*diag(ones(1,126)) ;
150 f = [zeros(1,126) C_lp_smoothness*ones(1,126) ones(1,N)/N] ; % <res, f'>
151
152
153
154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155 % qp_solve: LB <= res <= UB
156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
157 INF = 1e20 ;
158 LB = [-INF*ones(1,126) zeros(1,126) zeros(1,N)] ; % lower bound for res
159 UB = [ INF*ones(1,2*126) INF*ones(1,N)] ; % upper bound for res
160
161 """