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.old_w = cb.matrix([0.0] ,(1,1))
  53
  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+'))
  58
  59    def plog(self,string):
  60       if self.logfh != None:
  61          self.logfh.write(string)
  62          self.logfh.flush()
  63       else:
  64          print string
  65
  66    def addConstraint(self, energy_deltas, idx):
  67       energy_deltas = cb.matrix(energy_deltas)
  68       loss = 1.0
  69       useMarginRescaling = True
  70       self.currentIdx = idx
  71       self.currentLoss = loss
  72       self.energy_deltas = cb.matrix(energy_deltas)
  73       scalar_prod = 0.0
  74
  75       energy_deltas = energy_deltas.T
  76
  77       assert useMarginRescaling
  78
  79       if loss < self.EPSILON:
  80           return False
  81
  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')
  89                  return False
  90          else:
  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')
  94                  return False
  95
  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,:]
  99
 100       self.newG[-1,0:self.numFeatures] = -1.0 * energy_deltas
 101
 102       if useMarginRescaling:
 103           # margin rescaling
 104           self.newG[-1,self.numFeatures+idx] = -1.0
 105       else:
 106           # slack rescaling
 107           self.newG[-1,self.numFeatures+idx] = -1.0 / loss
 108       self.G = self.newG
 109
 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]
 113
 114       if useMarginRescaling:
 115           # margin rescaling
 116           self.new_h[-1,0] = -1.0 * loss
 117       else:
 118           # slack rescaling
 119           self.new_h[-1,0] = -1.0
 120       self.h = self.new_h
 121
 122       self.numConstraints = self.G.size[0]
 123       self.plog('number of constraints is %d\n' % self.numConstraints)
 124
 125       assert self.G.size[0] == self.h.size[0]
 126
 127       return True
 128
 129    def solve(self):
 130       d = cs.qp(self.P,self.q,self.G,self.h)
 131       self.solver_runs += 1
 132
 133       self.plog('Return status of solver: ' + d['status'])
 134       new_w = d['x']
 135       self.old_w = new_w
 136       self.plog('slacks are:\n')
 137
 138       slacks_string = ''
 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'
 142
 143       self.plog(slacks_string+'\n')
 144       self.plog('end of slacks\n')
 145
 146       obj_value = 0.5 * new_w.T*self.P*new_w + self.q.T*new_w
 147       obj_value = obj_value[0]
 148
 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))
 154
 155       print 'SIQP_CVXOPT: Objective is: %e'%obj_value
 156       self.plog('SIQP_CVXOPT: Objective is: %e\n'%obj_value)
 157
 158       print self.old_objective_value
 159       print obj_value
 160       assert self.old_objective_value <= obj_value + self.EPSILON
 161
 162       self.old_objective_value = obj_value
 163
 164       #scalar_prod = self.currentLoss + (self.energy_deltas * new_w[0:self.numFeatures,0])
 165       #assert ( scalar_prod <= new_w[self.numFeatures+self.currentIdx] )
 166
 167       return obj_value,new_w[0:self.numFeatures,0],new_w[self.numFeatures:,0]