← Back to branch summary

~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

« back to all changes in this revision

Viewing changes to fieldlp.cpp

  • Committer: Package Import Robot
  • Author(s): Cédric Boutillier
  • Date: 2013-07-09 10:44:01 UTC
  • mfrom: (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130709104401-5q66ozz5j5af0dak
Tags: 0.5+dfsg-3
* Upload to unstable.
* modify remove_failing_tests_on_32bits.patch to replace command of
  0009RenderStairCase test with an empty one instead of deleting it.
* remove lintian override about spelling error

Show diffs side-by-side

added added

removed removed

Lines of Context:
fieldlp.cpp
 
1
#include "fieldlp.h"
 
2
 
 
3
FieldVector FieldLP::subvector(FieldVector const &v, set<int> const &b)
 
4
{
 
5
  FieldVector ret(v.getField(),b.size());
 
6
  int I=0;
 
7
  for(set<int>::const_iterator i=b.begin();i!=b.end();i++,I++)
 
8
    ret[I]=v[*i];
 
9
 
 
10
  return ret;
 
11
}
 
12
 
 
13
 
 
14
FieldMatrix FieldLP::submatrix(FieldMatrix const &m, set<int> const &b)
 
15
{
 
16
  FieldMatrix ret(m.getField(),b.size(),m.getWidth());
 
17
  int I=0;
 
18
  for(set<int>::const_iterator i=b.begin();i!=b.end();i++,I++)
 
19
    ret[I]=m[*i];
 
20
 
 
21
  return ret;
 
22
}
 
23
 
 
24
FieldVector FieldLP::edgeDirection(int i)const
 
25
{
 
26
  set<int> basis2=basis;
 
27
  assert(basis2.count(i));
 
28
  basis2.erase(i);
 
29
  FieldMatrix A2=submatrix(A,basis2);
 
30
  FieldMatrix ker=A2.reduceAndComputeKernel();
 
31
  assert(ker.getHeight()==1);
 
32
  FieldVector e=ker[0];
 
33
  if(dot(A[i],e).sign()<0)return e;
 
34
  return theField.zHomomorphism(-1)*e;
 
35
}
 
36
 
 
37
bool FieldLP::isImprovingDirection(int i)const
 
38
{
 
39
  return dot(edgeDirection(i),w).sign()>0;
 
40
}
 
41
 
 
42
FieldElement FieldLP::improvement(int i, int &newBasisMember)const
 
43
{
 
44
  FieldVector e=edgeDirection(i);
 
45
  FieldElement ew=dot(e,w);
 
46
  FieldElement ret(theField);
 
47
  bool first=true;
 
48
  newBasisMember=-1;
 
49
  for(int s=0;s<A.getHeight();s++)//This is _the_ line where an anticycling rule is implemented
 
50
  //  for(int s=A.getHeight()-1;s>=0;s--)//This is _the_ line where an anticycling rule is implemented
 
51
    {
 
52
      FieldElement denominator=dot(A[s],e);
 
53
      if(denominator.sign()>0)
 
54
        {
 
55
          FieldElement imp=(b[s]-dot(A[s],x))*denominator.inverse()*ew;
 
56
          if(first){ret=imp;newBasisMember=s;first=false;}
 
57
          if((imp-ret).sign()<=0){newBasisMember=s;ret=imp;}
 
58
        }
 
59
    }
 
60
  return ret;
 
61
}
 
62
 
 
63
int FieldLP::step()
 
64
{
 
65
  FieldElement impBest(theField);
 
66
  int bestIndex=-1;
 
67
  int bestNewElement=-1;
 
68
  for(set<int>::const_iterator i=basis.begin();i!=basis.end();i++)
 
69
    {
 
70
      //  AsciiPrinter P(Stderr);
 
71
      //  edgeDirection(*i).print(P);
 
72
    if(isImprovingDirection(*i))
 
73
      {
 
74
        //fprintf(Stderr,"IMPROVING DIRECTION\n");
 
75
        int newBasisElement;
 
76
        FieldElement imp=improvement(*i,newBasisElement);
 
77
        //P.printFieldElement(imp);
 
78
        //P.printNewLine();
 
79
        if(newBasisElement==-1)return -1; //UNBOUNDED
 
80
        if((imp-impBest).sign()>=0)//needed for perturbation
 
81
          {
 
82
            impBest=imp;
 
83
            bestIndex=*i;
 
84
            bestNewElement=newBasisElement;
 
85
          }
 
86
      }
 
87
    }
 
88
  if(bestIndex==-1)
 
89
    return 0; //OPTIMAL
 
90
 
 
91
  set<int> basis2=basis;
 
92
  //fprintf(Stderr,"REMOVING %i INSERTING %i\n",bestIndex,bestNewElement);
 
93
  basis2.erase(bestIndex);
 
94
  basis2.insert(bestNewElement);
 
95
  setCurrentBasis(basis2);
 
96
  return 1;
 
97
}
 
98
 
 
99
FieldLP::FieldLP(FieldMatrix const &A_, FieldVector const &b_):
 
100
  A(A_),
 
101
  b(b_),
 
102
  theField(A_.getField()),
 
103
  x(A_.getField(),0),
 
104
  w(A.getField(),A_.getWidth())
 
105
{
 
106
  //  if(A.getHeight()!=b.size())fprintf(Stderr,"%i %i",A.getHeight(),b.size());
 
107
  assert(A.getHeight()==b.size());
 
108
}
 
109
 
 
110
 
 
111
void FieldLP::setObjectiveFunction(FieldVector const &w_)
 
112
{
 
113
  w=w_;
 
114
  assert(A.getWidth()==w.size());
 
115
}
 
116
 
 
117
 
 
118
FieldVector FieldLP::basisToPoint(set<int> const &basis)const
 
119
{
 
120
  AsciiPrinter P(Stdout);
 
121
  //  FieldMatrix AT=A.transposed();
 
122
  //  AT.printMatrix(P);
 
123
  FieldMatrix A2(A.getField(),basis.size(),A.getWidth());
 
124
  FieldVector b2(A.getField(),basis.size());
 
125
  int I=0;
 
126
  for(set<int>::const_iterator i=basis.begin();i!=basis.end();i++,I++)
 
127
    {
 
128
      A2[I]=A[*i];
 
129
      b2[I]=b[*i];
 
130
    }
 
131
  //A2.printMatrix(P);
 
132
  //b2.print(P);
 
133
  FieldMatrix s=A2.solver();
 
134
  //fprintf(Stderr,"%i %i\n",s.getWidth(),b2.size()+A.getHeight());
 
135
  //s.printMatrix(P);
 
136
  //concatenation(b2,FieldVector(b2.getField(),A.getWidth())).print(P);
 
137
  FieldVector r=s.canonicalize(concatenation(b2,FieldVector(b2.getField(),A.getWidth())));
 
138
  //r.print(P);
 
139
  
 
140
  assert(r.subvector(0,basis.size()).isZero());
 
141
 
 
142
  return r.subvector(basis.size(),r.size());
 
143
}
 
144
 
 
145
 
 
146
void FieldLP::setCurrentBasis(set<int> const &basis_)
 
147
{
 
148
  basis=basis_;
 
149
  x=basisToPoint(basis);
 
150
}
 
151
 
 
152
 
 
153
void FieldLP::print(Printer &P)const
 
154
{
 
155
  P.printString("LPprogram\n");
 
156
  P.printString("A:\n");
 
157
  A.printMatrix(P);
 
158
  P.printString("w:\n");
 
159
  w.print(P);
 
160
  P.printString("b:\n");
 
161
  b.print(P);
 
162
  P.printString("current basis:\n");
 
163
  for(set<int>::const_iterator i=basis.begin();i!=basis.end();i++)P.printInteger(*i);  
 
164
  P.printString("x:\n");
 
165
  x.print(P);  
 
166
}
 
167
 
 
168
FieldMatrix FieldLP::computeLinealitySpace()
 
169
{
 
170
  FieldMatrix temp=A;
 
171
  return temp.reduceAndComputeKernel();
 
172
}
 
173
 
 
174
 
 
175
FieldLP FieldLP::buildLPForFeasibilityCheck()
 
176
{
 
177
  assert(computeLinealitySpace().getHeight()==0);
 
178
  set<int> S;
 
179
  {
 
180
    FieldMatrix temp=A.flipped().transposed();//important that inequalities which do not appear in basis set are
 
181
    temp.reduce();
 
182
    int i=-1;
 
183
    int j=-1;
 
184
    while(temp.nextPivot(i,j))
 
185
      {
 
186
        S.insert(A.getHeight()-1-j);//perturbed furthest
 
187
        //fprintf(Stderr,"INSERTING %i\n",A.getHeight()-1-j);
 
188
      }
 
189
  }
 
190
  FieldMatrix AS=submatrix(A,S);
 
191
  FieldMatrix ASSolver=AS.solver();
 
192
  FieldVector xs2=ASSolver.canonicalize(concatenation(subvector(b,S),FieldVector(b.getField(),A.getWidth()))).subvector(S.size(),S.size()+A.getWidth());
 
193
 
 
194
 
 
195
 
 
196
  set<int> inequalitiesWithSlack;
 
197
  for(int i=0;i<A.getHeight();i++)
 
198
    {
 
199
      if(!S.count(i))
 
200
        if((dot(A[i],xs2)-b[i]).sign()>0)
 
201
          {
 
202
            inequalitiesWithSlack.insert(i);
 
203
            //fprintf(Stderr,"adding slack variable\n");
 
204
          }
 
205
    }
 
206
  FieldMatrix newA=combineOnTop(
 
207
                                combineLeftRight(
 
208
                                                 A,
 
209
                                                 FieldMatrix(
 
210
                                                             theField,
 
211
                                                             A.getHeight(),
 
212
                                                             inequalitiesWithSlack.size()
 
213
                                                             )
 
214
                                                 ),
 
215
                                FieldMatrix(
 
216
                                            theField,
 
217
                                            inequalitiesWithSlack.size(),
 
218
                                            A.getWidth()+inequalitiesWithSlack.size()
 
219
                                            )
 
220
                                );
 
221
  FieldVector newW(theField,inequalitiesWithSlack.size()+A.getWidth());
 
222
  int I=0;
 
223
  for(set<int>::const_iterator i=inequalitiesWithSlack.begin();i!=inequalitiesWithSlack.end();i++,I++)
 
224
    {
 
225
      newA[*i][I+A.getWidth()]=theField.zHomomorphism(-1);
 
226
      //      int sign=dot(A[*i],xs2).sign();
 
227
      newA[I+A.getHeight()][I+A.getWidth()]=theField.zHomomorphism(-1);//(sign>0)?-1:1);
 
228
      newW[I+A.getWidth()]=theField.zHomomorphism(-1);//(sign>0)?-1:1);
 
229
    }
 
230
  set<int> newBasis;
 
231
 
 
232
  //  for(int i=0;i<A.getHeight();i++)newBasis.insert(i);
 
233
  for(set<int>::const_iterator i=S.begin();i!=S.end();i++)newBasis.insert(*i);
 
234
  for(set<int>::const_iterator i=inequalitiesWithSlack.begin();i!=inequalitiesWithSlack.end();i++)newBasis.insert(*i);
 
235
 
 
236
  FieldLP ret(newA,concatenation(b,FieldVector(theField,inequalitiesWithSlack.size())));
 
237
  ret.setObjectiveFunction(newW);
 
238
  ret.setCurrentBasis(newBasis);
 
239
 
 
240
  return ret;
 
241
}
 
242
 
 
243
 
 
244
bool FieldLP::findFeasibleBasis()
 
245
{
 
246
  FieldLP A2=buildLPForFeasibilityCheck();
 
247
 
 
248
      AsciiPrinter Q(Stdout);
 
249
  int status;
 
250
  do
 
251
    {
 
252
      //A2.print(Q);
 
253
      status=A2.step();
 
254
    }
 
255
  while(status==1);
 
256
  fprintf(Stderr,status?"LP is unbounded.\n":"Optimal solution found.\n");
 
257
  //A2.print(Q);
 
258
 
 
259
  set<int> newBasis;
 
260
  for(set<int>::const_iterator i=A2.basis.begin();i!=A2.basis.end();i++)
 
261
    {
 
262
      if(*i<A.getHeight())
 
263
        newBasis.insert(*i);
 
264
    }
 
265
 
 
266
  if(newBasis.size()!=A.getWidth())return false;
 
267
  setCurrentBasis(newBasis);
 
268
  return true;
 
269
}
 
270
 
 
271
 
 
272
FieldLP FieldLP::withNoLineality()
 
273
{
 
274
  set<int> S;
 
275
  {
 
276
    FieldMatrix temp=A;
 
277
    temp.reduce();
 
278
    int i=-1;
 
279
    int j=-1;
 
280
    while(temp.nextPivot(i,j))S.insert(j);
 
281
  }
 
282
  //  FieldMatrix(theField,A.getHeight(),S.size());
 
283
  FieldLP ret(submatrix(A.transposed(),S).transposed(),b);
 
284
 
 
285
  //fprintf(Stderr,"New height=%i New width=%i\n",ret.A.getHeight(),ret.A.getWidth());
 
286
 
 
287
  //  AsciiPrinter P(Stderr);
 
288
  //ret.print(P);
 
289
 
 
290
  if(w.size())ret.setObjectiveFunction(subvector(w,S));
 
291
 
 
292
  return ret;
 
293
}