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.
6 # Written (W) 2008 Fabio De Bona
7 # Copyright (C) 2008 Max-Planck-Society
16 import cvxopt
.base
as cb
17 import cvxopt
.solvers
as cs
24 cs
.options
['show_progress'] = True # default: True
26 cs
.options
['show_progress'] = False # default: True
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
36 class SIQPSolver(SIQP
):
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.
44 def __init__(self
,fSize
,numExamples
,c
,logfh
,settings
):
45 SIQP
.__init
__(self
,fSize
,numExamples
,c
,settings
)
47 self
.numConstraints
= 0
49 self
.old_objective_value
= -(sys
.maxint
-1)
52 self
.old_w
= cb
.matrix([0.0] ,(1,1))
54 #self.P.tofile(open('matrix_P_%d_%d'%(self.P.size[0],self.P.size[1]),'w+'))
55 #self.q.tofile(open('matrix_q_%d_%d'%(self.q.size[0],self.q.size[1]),'w+'))
56 #self.G.tofile(open('matrix_G_%d_%d'%(self.G.size[0],self.G.size[1]),'w+'))
57 #self.h.tofile(open('matrix_h_%d_%d'%(self.h.size[0],self.h.size[1]),'w+'))
59 def plog(self
,string
):
60 if self
.logfh
!= None:
61 self
.logfh
.write(string
)
66 def addConstraint(self
, energy_deltas
, idx
):
67 energy_deltas
= cb
.matrix(energy_deltas
)
69 useMarginRescaling
= True
71 self
.currentLoss
= loss
72 self
.energy_deltas
= cb
.matrix(energy_deltas
)
75 energy_deltas
= energy_deltas
.T
77 assert useMarginRescaling
79 if loss
< self
.EPSILON
:
82 if self
.old_w
.size
!= (1,1):
83 scalar_prod
= energy_deltas
* self
.old_w
[0:self
.numFeatures
,0]
84 old_slack
= self
.old_w
[self
.numFeatures
+idx
,0]
85 if useMarginRescaling
:
86 if scalar_prod
[0] + old_slack
> loss
: # leave some room for inaccuracies
87 print 'Current example already fulfills constraint.'
88 self
.plog('Current example already fulfills constraint.\n')
91 self
.plog("constraint at current solution: %f <= -1+%f/%f = %f" %(scalar_prod
[0], old_slack
, loss
, -1 + old_slack
/loss
) )
92 if scalar_prod
[0] < -1+old_slack
/loss
+ 1e-6: # leave some room for inaccuracies
93 self
.plog('Current example already fulfills constraint.\n')
96 rows
,cols
= self
.G
.size
97 self
.newG
= cb
.matrix(0.0,(rows
+1,cols
))
98 self
.newG
[0:rows
,:] = self
.G
[0:rows
,:]
100 self
.newG
[-1,0:self
.numFeatures
] = -1.0 * energy_deltas
102 if useMarginRescaling
:
104 self
.newG
[-1,self
.numFeatures
+idx
] = -1.0
107 self
.newG
[-1,self
.numFeatures
+idx
] = -1.0 / loss
110 rows
,cols
= self
.h
.size
111 self
.new_h
= cb
.matrix(0.0,(rows
+1,cols
))
112 self
.new_h
[0:rows
,0] = self
.h
[0:rows
,0]
114 if useMarginRescaling
:
116 self
.new_h
[-1,0] = -1.0 * loss
119 self
.new_h
[-1,0] = -1.0
122 self
.numConstraints
= self
.G
.size
[0]
123 self
.plog('number of constraints is %d\n' % self
.numConstraints
)
125 assert self
.G
.size
[0] == self
.h
.size
[0]
130 d
= cs
.qp(self
.P
,self
.q
,self
.G
,self
.h
)
131 self
.solver_runs
+= 1
133 self
.plog('Return status of solver: ' + d
['status'])
136 self
.plog('slacks are:\n')
139 for i
in range(self
.numFeatures
,new_w
.size
[0]):
140 slacks_string
+= ('%e ' % new_w
[i
])
141 assert new_w
[i
] >= 0, 'Error found a negative slack variable'
143 self
.plog(slacks_string
+'\n')
144 self
.plog('end of slacks\n')
146 obj_value
= 0.5 * new_w
.T
*self
.P
*new_w
+ self
.q
.T
*new_w
147 obj_value
= obj_value
[0]
149 # evalute parts of the objective
150 regularizationPart
= (0.5 *new_w
.T
*self
.P
*new_w
)[0]
151 lossPart
= (self
.q
.T
*new_w
)[0]
152 assert lossPart
>= 0.0, 'Error we have a negative loss contribution'
153 self
.plog('Parts of objective: %e %e\n'%(regularizationPart
,lossPart
))
155 print 'SIQP_CVXOPT: Objective is: %e'%obj
_value
156 self
.plog('SIQP_CVXOPT: Objective is: %e\n'%obj
_value
)
158 print self
.old_objective_value
160 assert self
.old_objective_value
<= obj_value
+ self
.EPSILON
162 self
.old_objective_value
= obj_value
164 #scalar_prod = self.currentLoss + (self.energy_deltas * new_w[0:self.numFeatures,0])
165 #assert ( scalar_prod <= new_w[self.numFeatures+self.currentIdx] )
167 return obj_value
,new_w
[0:self
.numFeatures
,0],new_w
[self
.numFeatures
:,0]