qpalma/SIQP.py

   1 #!/usr/bin/env python
   2 # -*- coding:latin-1 -*-
   3
   4 import sys
   5 import os
   6 import os.path
   7 import logging
   8 import math
   9 import pdb
  10
  11 import cvxopt.base as cb
  12
  13 import logging
  14 logging.basicConfig(level=logging.DEBUG,format='%(levelname)s %(message)s')
  15
  16 import qpalma.Configuration as Configuration
  17
  18 class SIQP:
  19    """
  20    This class is a wrapper for the cvxopt/cplex qp solvers.
  21    It has the purpose to support an iterative addition of
  22    constraints to the original problem.
  23
  24    We want to solve a problem of the form
  25
  26    min  1/2 * C * x' * P * x + q' * x
  27     x
  28
  29     s.t.
  30       < w,f+ > - < w,f- > >= 1 - xi_i    for all i, for all f-
  31
  32    where x is a vector storing information on the parameters w_i and the slacks xi_i
  33    so x = [w;xi]
  34
  35    """
  36
  37    def __init__(self,fSize,numExamples,c):
  38       assert fSize > 0, 'You have to have at least one feature!'
  39       assert numExamples > 0, 'You have to have at least one example!'
  40       self.numFeatures = fSize
  41       self.numExamples = numExamples
  42       self.C = c
  43       self.EPSILON = 10e-15
  44
  45       logging.debug("Starting SIQP solver with %d examples. A feature space dim. of %d.", numExamples,fSize)
  46       logging.debug("Regularization parameters are C=%f."%(self.C))
  47
  48       self.old_w = None
  49
  50       self.dimP = self.numFeatures + self.numExamples
  51       self.createRegularizer()
  52
  53    def createUnitRegularizer(self):
  54       # set regularizer to zero
  55       self.P = cb.matrix(0.0,(self.dimP,self.dimP))
  56       for i in range(self.numFeatures):
  57          self.P[i,i] = 1.0
  58       # end of zeroing regularizer
  59
  60    def createSmoothnessRegularizer(self):
  61       # set regularizer to zero
  62       self.P = cb.matrix(0.0,(self.dimP,self.dimP))
  63       for i in range(self.numFeatures):
  64          self.P[i,i] = 1.0
  65
  66       lengthSP    = Configuration.numLengthSuppPoints
  67       donSP       = Configuration.numDonSuppPoints
  68       accSP       = Configuration.numAccSuppPoints
  69       mmatrixSP   = Configuration.sizeMatchmatrix[0]\
  70       *Configuration.sizeMatchmatrix[1]
  71       numq        = Configuration.numQualSuppPoints
  72       totalQualSP = Configuration.totalQualSuppPoints
  73       numQPlifs = Configuration.numQualPlifs
  74
  75       #lengthGroupParam  = 5*1e-9
  76       #spliceGroupParam  = 1e-8
  77
  78       lengthGroupParam  = 5*1e-9
  79       spliceGroupParam  = 1e-9
  80
  81       self.P[0:lengthSP,0:lengthSP] *= lengthGroupParam
  82
  83       beg = lengthSP
  84       end = lengthSP+donSP+accSP
  85       self.P[beg:end,beg:end]       *= spliceGroupParam
  86
  87       regC = self.numFeatures / (1.0*self.numExamples)
  88
  89       #cfactor = 1e-8
  90       cfactor = 5e-9
  91
  92       for j in range(lengthSP-1):
  93          self.P[j,j+1]    = -cfactor
  94          self.P[j+1,j]    = -cfactor
  95          self.P[j,j]     += cfactor
  96          self.P[j+1,j+1] += cfactor
  97
  98       cfactor = 1e-9
  99       #cfactor = 5e-9
 100
 101       for j in range(lengthSP,lengthSP+donSP-1):
 102          cfactor = spliceGroupParam
 103          self.P[j,j+1] = -cfactor
 104          self.P[j+1,j] = -cfactor
 105          self.P[j,j] += cfactor
 106          self.P[j+1,j+1] += cfactor
 107
 108       for j in range(lengthSP+donSP,lengthSP+donSP+accSP-1):
 109          cfactor = spliceGroupParam
 110          self.P[j,j+1] = -cfactor
 111          self.P[j+1,j] = -cfactor
 112          self.P[j,j] += cfactor
 113          self.P[j+1,j+1] += cfactor
 114
 115       offset = lengthSP+donSP+accSP+mmatrixSP
 116       for k in range(numQPlifs):
 117          begin = offset+(k*numq)
 118          end   = offset+((k+1)*numq)
 119
 120          for j in range(begin,end-1):
 121             self.P[j,j+1] = -1.0
 122             self.P[j+1,j] = -1.0
 123             self.P[j,j] += 1.0
 124
 125       # 0.25 for each was already good
 126
 127       matchGroupParam   = 1.0
 128       qualityGroupParam = 1.0
 129
 130       # lengthGroupParam  = 0.005
 131       # spliceGroupParam  = 0.005
 132       # matchGroupParam   = 0.495
 133       # qualityGroupParam = 0.495
 134
 135       #lengthGroupParam  = lengthSP / (1.0*lengthSP+donSP+accSP)
 136       #spliceGroupParam  = donSP+accSP / (1.0*lengthSP+donSP+accSP)
 137       #matchGroupParam   = mmatrixSP/(1.0*mmatrixSP+totalQualSP)
 138       #qualityGroupParam = totalQualSP/(1.0*mmatrixSP+totalQualSP)
 139
 140
 141       #denom = (1.0*lengthSP+donSP+accSP+mmatrixSP+totalQualSP)
 142       #ratio1 = lengthSP+donSP+accSP / denom
 143       #ratio2 = mmatrixSP+totalQualSP / denom
 144
 145       #lengthGroupParam  *= ratio1
 146       #spliceGroupParam  *= ratio1
 147
 148       #matchGroupParam   *= ratio2
 149       #qualityGroupParam *= ratio2
 150
 151
 152       beg = lengthSP+donSP+accSP
 153       end = lengthSP+donSP+accSP+mmatrixSP
 154       self.P[beg:end,beg:end]       *= matchGroupParam
 155
 156       beg = lengthSP+donSP+accSP+mmatrixSP
 157       self.P[beg:-self.numExamples,beg:-self.numExamples]   *= qualityGroupParam
 158
 159       self.P[0:-self.numExamples,0:-self.numExamples]       *= regC*0.1
 160
 161
 162    def createRegularizer(self):
 163       #self.createUnitRegularizer()
 164       self.createSmoothnessRegularizer()
 165
 166       q_array = [0.0]*self.numFeatures + [1.0]*self.numExamples
 167       self.q = cb.matrix(q_array,(self.numFeatures+self.numExamples,1),'d')
 168       self.q = self.C * self.q
 169
 170       self.G = cb.matrix(0.0,(self.numExamples,self.numFeatures+self.numExamples))
 171       for i in range(self.numExamples):
 172          self.G[i,self.numFeatures+i] = -1
 173
 174       self.h = cb.matrix(0,(self.numExamples,1),'d')
 175
 176 if __name__ == '__main__':
 177    sq = SIQP(3,2,100.0)
 178
 179    """
 180    function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
 181    % function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
 182
 183    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 184    % qp_solve: min ( <res, f'> + 1/2 res' Q res)
 185    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 186    C_qp_penalties  = C ;    %was 0.01
 187    C_qp_smoothness = C ;    %was 0.01
 188    C_lp_smoothness = 0.00 ; %linear problem parameters
 189    C_qp_smallness  = 1e-6 ; %was 1e-6
 190    column_eps      = 1e-3 ;
 191
 192    Q=zeros(2*126+N) ; % 1/2 * res' * Q * res
 193    Q(1:90,1:90)           = C_qp_smallness*diag(ones(1,90)) ;
 194    Q(91:126,91:126)       = C_qp_penalties*diag(ones(1,36)) ;
 195    Q(127:2*126,127:2*126) = C_qp_smoothness*diag(ones(1,126)) ;
 196    f = [zeros(1,126) C_lp_smoothness*ones(1,126) ones(1,N)/N] ; % <res, f'>
 197
 198    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 199    % qp_solve: LB <= res <= UB
 200    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 201    INF = 1e20 ;
 202    LB = [-INF*ones(1,126) zeros(1,126) zeros(1,N)] ; % lower bound for res
 203    UB = [ INF*ones(1,2*126) INF*ones(1,N)] ; % upper bound for res
 204
 205    """