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