python/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 def asymCoords(t):
  17    size_upstream = t[0]
  18    size_downstream = t[1]
  19    assert size_upstream >= 1 and size_upstream <= 39
  20    assert size_downstream >= 1 and size_downstream <= 39
  21
  22    if size_downstream > size_upstream:
  23       asym_small = size_upstream
  24       asym_big   = size_downstream
  25    else:
  26       asym_small = size_downstream
  27       asym_big   = size_upstream
  28
  29    asym_offset = 0
  30    for r in range(asym_small):
  31       asym_offset += (40-r)-r;
  32
  33    return (asym_offset+asym_big)-1
  34
  35 def calcTuples(i,j):
  36    return filter(lambda t: t[0]<=t[1] and t[0]>0 and t[0]<40 and t[1]>0 and t[1]<40,[(i+1,j),(i-1,j),(i,j+1),(i,j-1)])
  37
  38 class SIQP:
  39    """
  40    This class is a wrapper for the cvxopt qp solver.
  41    It has the purpose to support an iterative addition of
  42    constraints to the original problem.
  43
  44    We want to solve a problem of the form
  45
  46    min  1/2 * C * x' * P * x + q' * x
  47     x
  48
  49     s.t.
  50       < w,f+ > - < w,f- > >= 1 - xi_i    for all i, for all f-
  51
  52    where x is a vector storing information on the parameters w_i and the slacks xi_i
  53    so x = [w;xi]
  54
  55    """
  56
  57    def __init__(self,fSize,numExamples,c):
  58       assert fSize > 0, 'You have to have at least one feature!'
  59       assert numExamples > 0, 'You have to have at least one example!'
  60       self.numFeatures = fSize
  61       self.numExamples = numExamples
  62       self.C = c
  63       self.EPSILON = 10e-15
  64
  65       logging.debug("Starting SIQP solver with %d examples. A feature space dim. of %d.", numExamples,fSize)
  66       logging.debug("Regularization parameters are C=%f."%(self.C))
  67
  68       self.old_w = None
  69
  70       self.dimP = self.numFeatures + self.numExamples
  71       self.createRegularizer()
  72
  73    def createUnitRegularizer(self):
  74       # set regularizer to zero
  75       self.P = cb.matrix(0.0,(self.dimP,self.dimP))
  76       for i in range(self.numFeatures):
  77          self.P[i,i] = 1.0
  78       # end of zeroing regularizer
  79
  80    def createRegularizer(self):
  81       self.createUnitRegularizer()
  82
  83       q_array = [0.0]*self.numFeatures + [1.0]*self.numExamples
  84       self.q = cb.matrix(q_array,(self.numFeatures+self.numExamples,1),'d')
  85       self.q = self.C * self.q
  86
  87       self.G = cb.matrix(0.0,(self.numExamples,self.numFeatures+self.numExamples))
  88       for i in range(self.numExamples):
  89          self.G[i,self.numFeatures+i] = -1
  90
  91       self.h = cb.matrix(0,(self.numExamples,1),'d')
  92
  93 if __name__ == '__main__':
  94    sq = SIQP(3,2,100.0)
  95
  96    """
  97    function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
  98    % function [C_qp_penalties,C_qp_smoothness,C_lp_smoothness,C_qp_smallness,column_eps,Q,f,LB,UB] = paths_create_qp(N,C)
  99
 100    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 101    % qp_solve: min ( <res, f'> + 1/2 res' Q res)
 102    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 103    C_qp_penalties  = C ;    %was 0.01
 104    C_qp_smoothness = C ;    %was 0.01
 105    C_lp_smoothness = 0.00 ; %linear problem parameters
 106    C_qp_smallness  = 1e-6 ; %was 1e-6
 107    column_eps      = 1e-3 ;
 108
 109    Q=zeros(2*126+N) ; % 1/2 * res' * Q * res
 110    Q(1:90,1:90)           = C_qp_smallness*diag(ones(1,90)) ;
 111    Q(91:126,91:126)       = C_qp_penalties*diag(ones(1,36)) ;
 112    Q(127:2*126,127:2*126) = C_qp_smoothness*diag(ones(1,126)) ;
 113    f = [zeros(1,126) C_lp_smoothness*ones(1,126) ones(1,N)/N] ; % <res, f'>
 114
 115
 116
 117    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 118    % qp_solve: LB <= res <= UB
 119    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 120    INF = 1e20 ;
 121    LB = [-INF*ones(1,126) zeros(1,126) zeros(1,N)] ; % lower bound for res
 122    UB = [ INF*ones(1,2*126) INF*ones(1,N)] ; % upper bound for res
 123
 124    """