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