+ added small program for intron-position comparison
[qpalma.git] / qpalma / SIQP_CVXOPT.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 import cvxopt.solvers as cs
13
14 import logging
15 logging.basicConfig(level=logging.DEBUG,format='%(levelname)s %(message)s')
16
17 from SIQP import SIQP
18
19 acc = 1e-7
20
21 if __debug__:
22 cs.options['show_progress'] = True # default: True
23 else:
24 cs.options['show_progress'] = False # default: True
25
26 cs.options['maxiters'] = 200 # default: 100
27 cs.options['abstol'] = acc # default: 1e-7
28 cs.options['reltol'] = acc # default: 1e-7
29 cs.options['feastol'] = acc # default: 1e-7
30 #cs.options['refinement'] = True # default: True
31
32 class SIQPSolver(SIQP):
33 """
34 This class is a wrapper for the cvxopt qp solver.
35 It has the purpose to support an iterative addition of
36 constraints to the original problem.
37
38 """
39
40 def __init__(self,fSize,numExamples,c,proto,run):
41 SIQP.__init__(self,fSize,numExamples,c,run)
42 self.numConstraints = 0
43 self.solver_runs = 0
44 self.old_objective_value = -(sys.maxint-1)
45 self.protocol = proto
46
47 #self.P.tofile(open('matrix_P_%d_%d'%(self.P.size[0],self.P.size[1]),'w+'))
48 #self.q.tofile(open('matrix_q_%d_%d'%(self.q.size[0],self.q.size[1]),'w+'))
49 #self.G.tofile(open('matrix_G_%d_%d'%(self.G.size[0],self.G.size[1]),'w+'))
50 #self.h.tofile(open('matrix_h_%d_%d'%(self.h.size[0],self.h.size[1]),'w+'))
51
52 def addConstraint(self, energy_deltas, idx):
53 energy_deltas = cb.matrix(energy_deltas)
54 loss = 1.0
55 useMarginRescaling = True
56 self.currentIdx = idx
57 self.currentLoss = loss
58 self.energy_deltas = cb.matrix(energy_deltas)
59 scalar_prod = 0.0
60
61 energy_deltas = energy_deltas.T
62
63 assert useMarginRescaling
64
65 if loss < self.EPSILON:
66 return False
67
68 if self.old_w != None:
69 scalar_prod = energy_deltas * self.old_w[0:self.numFeatures,0]
70 old_slack = self.old_w[self.numFeatures+idx,0]
71 if useMarginRescaling:
72 if scalar_prod[0] + old_slack > loss: # leave some room for inaccuracies
73 print 'Current example already fulfills constraint.'
74 print >> self.protocol, 'Current example already fulfills constraint.'
75 return False
76 else:
77 print >> self.protocol, "constraint at current solution: %f <= -1+%f/%f = %f" %(scalar_prod[0], old_slack, loss, -1 + old_slack/loss)
78 if scalar_prod[0] < -1+old_slack/loss + 1e-6: # leave some room for inaccuracies
79 print >> self.protocol, 'Current example already fulfills constraint.'
80 logging.debug('Current example already fulfills constraint.')
81 return False
82
83 rows,cols = self.G.size
84 self.newG = cb.matrix(0.0,(rows+1,cols))
85 self.newG[0:rows,:] = self.G[0:rows,:]
86
87 self.newG[-1,0:self.numFeatures] = -1.0 * energy_deltas
88
89 if useMarginRescaling:
90 # margin rescaling
91 self.newG[-1,self.numFeatures+idx] = -1.0
92 else:
93 # slack rescaling
94 self.newG[-1,self.numFeatures+idx] = -1.0 / loss
95 self.G = self.newG
96
97 rows,cols = self.h.size
98 self.new_h = cb.matrix(0.0,(rows+1,cols))
99 self.new_h[0:rows,0] = self.h[0:rows,0]
100
101 if useMarginRescaling:
102 # margin rescaling
103 self.new_h[-1,0] = -1.0 * loss
104 else:
105 # slack rescaling
106 self.new_h[-1,0] = -1.0
107 self.h = self.new_h
108
109 self.numConstraints = self.G.size[0]
110 print >> self.protocol, 'number of constraints is %d' % self.numConstraints
111
112 assert self.G.size[0] == self.h.size[0]
113
114 return True
115
116 def solve(self):
117 d = cs.qp(self.P,self.q,self.G,self.h)
118 self.solver_runs += 1
119
120 print >> self.protocol, 'Return status of solver: ' + d['status']
121 new_w = d['x']
122 self.old_w = new_w
123 print >> self.protocol, 'slacks are'
124 for i in range(self.numFeatures,new_w.size[0]):
125 print >> self.protocol, "%e " % new_w[i],
126 assert new_w[i] >= 0, 'Error found a negative slack variable'
127 print >> self.protocol, "end of slacks"
128
129 obj_value = 0.5 * new_w.T*self.P*new_w + self.q.T*new_w
130 obj_value = obj_value[0]
131
132 # evalute parts of the objective
133 regularizationPart = (0.5 *new_w.T*self.P*new_w)[0]
134 lossPart = (self.q.T*new_w)[0]
135 assert lossPart >= 0.0, 'Error we have a negative loss contribution'
136 print >> self.protocol, 'Parts of objective: %e %e'%(regularizationPart,lossPart)
137
138 print 'SIQP_CVXOPT: Objective is: %e'%obj_value
139 print >> self.protocol, 'SIQP_CVXOPT: Objective is: %e'%obj_value
140
141 print self.old_objective_value
142 print obj_value
143 assert self.old_objective_value <= obj_value + self.EPSILON
144
145 self.old_objective_value = obj_value
146
147 #scalar_prod = self.currentLoss + (self.energy_deltas * new_w[0:self.numFeatures,0])
148 #assert ( scalar_prod <= new_w[self.numFeatures+self.currentIdx] )
149
150 return obj_value,new_w[0:self.numFeatures,0],new_w[self.numFeatures:,0]
151
152 def test():
153 fh = open('siqp_cvxopt.log','w+')
154 numFeatures = 3
155 numExamples = 2
156 s = SIQPSolver(numFeatures,numExamples,100.0,fh)
157
158 d1 = cb.matrix([1,0,2],(numFeatures,1))
159 d2 = cb.matrix([1,1,0],(numFeatures,1))
160
161 s.addConstraint(d1,10.0,0,True)
162 s.addConstraint(d2,105.0,1,True)
163
164 w,slacks = s.solve()
165 print w
166 print slacks
167
168 print w.T * d1 + slacks[0]
169 print w.T * d2 + slacks[1]
170
171 del s
172 fh.close()
173
174 if __name__ == '__main__':
175 test()