+ added solver wrapper for cvxopt
authorfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Fri, 18 Apr 2008 09:12:12 +0000 (09:12 +0000)
committerfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Fri, 18 Apr 2008 09:12:12 +0000 (09:12 +0000)
+ modified parser

git-svn-id: http://svn.tuebingen.mpg.de/ag-raetsch/projects/QPalma@8659 e1793c9e-67f9-0310-80fc-b846ff1f7b36

qpalma/SIQP_CVXOPT.py [new file with mode: 0644]
qpalma/parsers.py

diff --git a/qpalma/SIQP_CVXOPT.py b/qpalma/SIQP_CVXOPT.py
new file mode 100644 (file)
index 0000000..c09340b
--- /dev/null
@@ -0,0 +1,175 @@
+#!/usr/bin/env python
+# -*- coding:latin-1 -*-
+
+import sys
+import os
+import os.path
+import logging
+import math
+import pdb
+
+import cvxopt.base as cb
+import cvxopt.solvers as cs
+
+import logging
+logging.basicConfig(level=logging.DEBUG,format='%(levelname)s %(message)s')
+
+from SIQP import SIQP
+
+acc = 1e-7
+
+if __debug__:
+   cs.options['show_progress'] = True   # default: True 
+else:
+   cs.options['show_progress'] = False  # default: True 
+
+cs.options['maxiters']      = 200      # default: 100
+cs.options['abstol']        = acc      # default: 1e-7
+cs.options['reltol']        = acc      # default: 1e-7
+cs.options['feastol']       = acc      # default: 1e-7
+#cs.options['refinement']    = True     # default: True 
+
+class SIQPSolver(SIQP):
+   """
+   This class is a wrapper for the cvxopt qp solver.
+   It has the purpose to support an iterative addition of
+   constraints to the original problem.
+   
+   """
+
+   def __init__(self,fSize,numExamples,c,proto,run):
+      SIQP.__init__(self,fSize,numExamples,c,run)
+      self.numConstraints = 0
+      self.solver_runs = 0
+      self.old_objective_value = -(sys.maxint-1)
+      self.protocol = proto
+
+      #self.P.tofile(open('matrix_P_%d_%d'%(self.P.size[0],self.P.size[1]),'w+'))
+      #self.q.tofile(open('matrix_q_%d_%d'%(self.q.size[0],self.q.size[1]),'w+'))
+      #self.G.tofile(open('matrix_G_%d_%d'%(self.G.size[0],self.G.size[1]),'w+'))
+      #self.h.tofile(open('matrix_h_%d_%d'%(self.h.size[0],self.h.size[1]),'w+'))
+
+   def addConstraint(self, energy_deltas, idx):
+      energy_deltas = cb.matrix(energy_deltas)
+      loss = 1.0
+      useMarginRescaling = True
+      self.currentIdx = idx
+      self.currentLoss = loss
+      self.energy_deltas = cb.matrix(energy_deltas)
+      scalar_prod = 0.0
+
+      energy_deltas = energy_deltas.T
+
+      assert useMarginRescaling
+
+      if loss < self.EPSILON:
+          return False
+      
+      if self.old_w != None:
+         scalar_prod = energy_deltas * self.old_w[0:self.numFeatures,0]
+         old_slack = self.old_w[self.numFeatures+idx,0]
+         if useMarginRescaling:
+             if scalar_prod[0] + old_slack > loss: # leave some room for inaccuracies
+                 print 'Current example already fulfills constraint.'
+                 print >> self.protocol,  'Current example already fulfills constraint.'
+                 return False
+         else:
+             print >> self.protocol, "constraint at current solution: %f <= -1+%f/%f = %f" %(scalar_prod[0], old_slack, loss, -1 + old_slack/loss) 
+             if scalar_prod[0] < -1+old_slack/loss + 1e-6: # leave some room for inaccuracies
+                 print >> self.protocol,  'Current example already fulfills constraint.'
+                 logging.debug('Current example already fulfills constraint.')
+                 return False
+
+      rows,cols = self.G.size
+      self.newG = cb.matrix(0.0,(rows+1,cols))
+      self.newG[0:rows,:] = self.G[0:rows,:]
+
+      self.newG[-1,0:self.numFeatures] = -1.0 * energy_deltas
+
+      if useMarginRescaling:
+          # margin rescaling
+          self.newG[-1,self.numFeatures+idx] = -1.0
+      else:
+          # slack rescaling
+          self.newG[-1,self.numFeatures+idx] = -1.0 / loss
+      self.G = self.newG
+         
+      rows,cols = self.h.size
+      self.new_h = cb.matrix(0.0,(rows+1,cols))
+      self.new_h[0:rows,0] = self.h[0:rows,0]
+
+      if useMarginRescaling:
+          # margin rescaling
+          self.new_h[-1,0] = -1.0 * loss
+      else:
+          # slack rescaling
+          self.new_h[-1,0] = -1.0
+      self.h = self.new_h
+
+      self.numConstraints = self.G.size[0]
+      print >> self.protocol, 'number of constraints is %d' % self.numConstraints
+
+      assert self.G.size[0] == self.h.size[0]
+
+      return True
+
+   def solve(self):
+      d = cs.qp(self.P,self.q,self.G,self.h)
+      self.solver_runs += 1
+
+      print >> self.protocol, 'Return status of solver: ' + d['status']
+      new_w = d['x']
+      self.old_w = new_w
+      print >> self.protocol, 'slacks are'
+      for i in range(self.numFeatures,new_w.size[0]):
+         print >> self.protocol, "%e " % new_w[i],
+         assert new_w[i] >= 0, 'Error found a negative slack variable'
+      print >> self.protocol, "end of slacks"
+
+      obj_value = 0.5 * new_w.T*self.P*new_w + self.q.T*new_w
+      obj_value = obj_value[0]
+      
+      # evalute parts of the objective
+      regularizationPart = (0.5 *new_w.T*self.P*new_w)[0]
+      lossPart = (self.q.T*new_w)[0]
+      assert lossPart >= 0.0, 'Error we have a negative loss contribution'
+      print >> self.protocol, 'Parts of objective: %e %e'%(regularizationPart,lossPart)
+
+      print 'SIQP_CVXOPT: Objective is: %e'%obj_value
+      print >> self.protocol, 'SIQP_CVXOPT: Objective is: %e'%obj_value
+
+      print self.old_objective_value 
+      print obj_value 
+      assert self.old_objective_value <= obj_value + self.EPSILON
+
+      self.old_objective_value = obj_value
+
+      #scalar_prod = self.currentLoss + (self.energy_deltas * new_w[0:self.numFeatures,0])
+      #assert ( scalar_prod <= new_w[self.numFeatures+self.currentIdx] )
+
+      return obj_value,new_w[0:self.numFeatures,0],new_w[self.numFeatures:,0]
+
+def test():
+   fh = open('siqp_cvxopt.log','w+')
+   numFeatures = 3
+   numExamples = 2
+   s = SIQPSolver(numFeatures,numExamples,100.0,fh)
+
+   d1 = cb.matrix([1,0,2],(numFeatures,1))
+   d2 = cb.matrix([1,1,0],(numFeatures,1))
+
+   s.addConstraint(d1,10.0,0,True)
+   s.addConstraint(d2,105.0,1,True)
+
+   w,slacks = s.solve()
+   print w
+   print slacks
+
+   print w.T * d1 + slacks[0]
+   print w.T * d2 + slacks[1]
+
+   del s
+   fh.close()
+
+if __name__ == '__main__':
+   test()
index c22a75f..62bc029 100644 (file)
@@ -41,7 +41,12 @@ class FilteredReadParser(ReadParser):
       read_nr,chr,strand,seq,splitpos,read_size,prb,cal_prb,chastity,gene_id,p_start,exon_stop,exon_start,p_stop
 
       """
-      id,chr,strand,seq,splitpos,read_size,prb,cal_prb,chastity,gene_id,p_start,exon_stop,exon_start,p_stop,true_cut = line.split()
+      try:
+         id,chr,strand,seq,splitpos,read_size,prb,cal_prb,chastity,gene_id,p_start,exon_stop,exon_start,p_stop,true_cut = line.split()
+      except:
+         id,chr,strand,seq,splitpos,read_size,prb,cal_prb,chastity,gene_id,p_start,exon_stop,exon_start,p_stop = line.split()
+         true_cut = -1
+
       splitpos = int(splitpos)
       read_size = int(read_size)