--- /dev/null
+#!/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()
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)
projects / qpalma.git / commitdiff