1'''
2As a part of ``cylp.python.pivots`` it implements the positive edge
3pivot selection rule.
4'''
5
6from __future__ import print_function
7import random
8import numpy as np
9from cylp.cy import CyCoinIndexedVector
10from cylp.cy.CyClpSimplex import cydot
11from .PivotPythonBase import PivotPythonBase
12
13
14class PositiveEdgePivot(PivotPythonBase):
15    '''
16    Positive Edge pivot rule implementation.
17
18    .. _custom-pivot-usage:
19
20    **Usage**
21
22    >>> from cylp.cy import CyClpSimplex
23    >>> from cylp.py.pivots import PositiveEdgePivot
24    >>> from cylp.py.pivots.PositiveEdgePivot import getMpsExample
25    >>> # Get the path to a sample mps file
26    >>> f = getMpsExample()
27    >>> s = CyClpSimplex()
28    >>> s.readMps(f)  # Returns 0 if OK
29    0
30    >>> pivot = PositiveEdgePivot(s)
31    >>> s.setPivotMethod(pivot)
32    >>> s.primal()
33    'optimal'
34    >>> round(s.objectiveValue, 5)
35    2520.57174
36
37    '''
38
39    def __init__(self, clpModel, EPSILON=10 ** (-7)):
40        self.clpModel = clpModel
41        self.dim = self.clpModel.nRows + self.clpModel.nCols
42
43        self.isDegenerate = False
44
45        # Create some numpy arrays here ONCE to prevent memory
46        # allocation at each iteration
47        self.aColumn = CyCoinIndexedVector()
48        self.aColumn.reserve(self.dim)
49        self.w = CyCoinIndexedVector()
50        self.w.reserve(self.clpModel.nRows)
51
52        self.rhs = np.empty(self.clpModel.nRows, dtype=np.double)
53        self.EPSILON = EPSILON
54        self.lastUpdateIteration = 0
55
56    def updateP(self):
57        '''Finds constraints with abs(rhs) <=  epsilon and put
58        their indices in "z"
59        '''
60        s = self.clpModel
61        nRows = s.nRows
62
63        rhs = self.rhs
64        s.getRightHandSide(rhs)
65
66        #self.p = np.where(np.abs(rhs) > self.EPSILON)[0]
67        self.z = np.where(np.abs(rhs) <= self.EPSILON)[0]
68
69        #print('degeneracy level : ', (len(self.z)) / float(nRows)))
70        self.isDegenerate = (len(self.z) > 0)
71
72    def updateW(self):
73        '''Sets "w" to be a vector of random vars with "0"
74        at indices defined in "p"
75        Note that vectorTimesB_1 changes "w"
76        '''
77        self.updateP()
78        self.w.clear()
79        self.w[self.z] = np.random.random(len(self.z))
80        s = self.clpModel
81        s.vectorTimesB_1(self.w)
82
83        self.lastUpdateIteration = s.iteration
84
85    def random(self):
86        'Defines how random vector "w" components are generated'
87        return random.random()
88
89    def isCompatible(self, varInd):
90        if not self.isDegenerate:
91            return False
92        s = self.clpModel
93        s.getACol(varInd, self.aColumn)
94
95        return abs(cydot(self.aColumn, self.w)) < self.EPSILON
96
97    def checkVar(self, i):
98        return self.isCompatible(i)
99
100    def pivotColumn(self, updates, spareRow1, spareRow2, spareCol1, spareCol2):
101        self.updateReducedCosts(updates, spareRow1, spareRow2, spareCol1, spareCol2)
102        s = self.clpModel
103        rc = s.reducedCosts
104
105        tol = s.dualTolerance
106        indicesToConsider = np.where(s.varNotFlagged & s.varNotFixed &
107                                     s.varNotBasic &
108                                     (((rc > tol) & s.varIsAtUpperBound) |
109                                     ((rc < -tol) & s.varIsAtLowerBound) |
110                                     s.varIsFree))[0]
111
112        rc2 = abs(rc[indicesToConsider])
113
114        maxRc = maxCompRc = maxInd = maxCompInd = -1
115
116        if self.isDegenerate:
117            w = self.w.elements
118            compatibility = np.zeros(s.nCols + s.nRows, dtype=np.double)
119            if len(indicesToConsider) > 0:
120                s.transposeTimesSubsetAll(indicesToConsider.astype(np.int64),
121                                          w, compatibility)
122            comp_varInds = indicesToConsider[np.where(abs(
123                                    compatibility[indicesToConsider]) <
124                                    self.EPSILON)[0]]
125
126            comp_rc = abs(rc[comp_varInds])
127            if len(comp_rc) > 0:
128                maxCompInd = comp_varInds[np.argmax(comp_rc)]
129                maxCompRc = rc[maxCompInd]
130
131        if len(rc2) > 0:
132            maxInd = indicesToConsider[np.argmax(rc2)]
133            maxRc = rc[maxInd]
134
135        del rc2
136        if maxCompInd != -1 and abs(maxCompRc) > 0.1 * abs(maxRc):
137            return maxCompInd
138        self.updateW()
139        return maxInd
140
141    def saveWeights(self, model, mode):
142        self.clpModel = model
143
144    def isPivotAcceptable(self):
145        return True
146
147
148def getMpsExample():
149    import os
150    import inspect
151    cylpDir = os.environ['CYLP_SOURCE_DIR']
152    return os.path.join(cylpDir, 'cylp', 'input', 'p0033.mps')
153
154