00001
00002
00003
00004 #include "matrices.h"
00005 #include "gcd.h"
00006
00007 #ifndef MOTZKIN_H_
00008 #define MOTZKIN_H_
00009
00014
00018
00022
00024 const int mm_no_modification = 0;
00025 const int mm_min_modification = 1;
00026 const int mm_max_modification = 2;
00027 const int mm_no_movie = 0;
00028 const int mm_movie = 1;
00029
00030
00035
00036 #ifdef no_default_arguments
00037 template <class int_item>
00038 void skeleton(matrix<int_item> & a, matrix<int_item> & f,
00039 matrix<int_item> & q, matrix<int_item> & e,
00040 int modification,
00041 int movie,
00042 ostr& movie_stream);
00043 #else
00044 template <class int_item>
00045 void skeleton(matrix<int_item> & a, matrix<int_item> & f,
00046 matrix<int_item> & q, matrix<int_item> & e,
00047 int modification = mm_no_modification,
00048 int movie = mm_no_movie,
00049 ostr& movie_stream = cout);
00050 #endif
00051
00052
00053
00054
00055
00056 template <class int_item>
00057 class motzkin_burger_algorithm
00058 {
00059 public:
00060 motzkin_burger_algorithm(matrix<int_item>&a, matrix<int_item>&f,
00061 matrix<int_item>&q, matrix<int_item>&e,
00062 int modification, int movie, ostr& movie_stream)
00063 : a(a), f(f), q(q), e(e),
00064 modification(modification), movie(movie), movie_stream(movie_stream){};
00065 void run();
00066
00067 private:
00068 matrix<int_item> & a;
00069 matrix<int_item> & f;
00070 matrix<int_item> & q;
00071 matrix<int_item> & e;
00072 int modification;
00073 int movie;
00074 ostr & movie_stream;
00075
00076 void gauss();
00077
00078 int reroof(index row);
00079 int balance(index j_plus, index j_minus);
00080 index min_modification();
00081 index max_modification();
00082 index fun(index j);
00083
00084 void movie_print3();
00085
00086 vector<index>common_zero;
00087 index n, m, r, current_ost, current_column;
00088 };
00089
00090 #ifdef no_default_arguments
00091 template <class int_item>
00092 void skeleton(matrix<int_item> & a, matrix<int_item> & f,
00093 matrix<int_item> & q, matrix<int_item> & e,
00094 int modification,
00095 int movie,
00096 ostr& movie_stream)
00097 #else
00098 template <class int_item>
00099 void skeleton(matrix<int_item> & a, matrix<int_item> & f,
00100 matrix<int_item> & q, matrix<int_item> & e,
00101 int modification,
00102 int movie,
00103 ostr& movie_stream)
00104 #endif
00105 {
00106 motzkin_burger_algorithm<int_item>
00107 alg(a, f, q, e, modification, movie, movie_stream);
00108 alg.run();
00109 }
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121 template <class int_item>
00122 void motzkin_burger_algorithm<int_item>::gauss()
00123 {
00124 index i, j;
00125 index m = a.get_m();
00126 index n = a.get_n();
00127
00128 f.resize(n, n);
00129 e.resize(0, n);
00130
00131
00132
00133 for (i = 0; i < n; i++)
00134 for (j = 0; j < n; j++)
00135 f(i, j) = 0;
00136 for (i = 0; i < n; i++)
00137 f(i, i) = 1;
00138
00139 q = transpose(a);
00140
00141 if (movie)
00142 {
00143 movie_stream << "1. Gaussian algorithm\n"
00144 << " ******************\n";
00145 movie_print3();
00146 }
00147
00148
00149 i = 0;
00150 while (i < q.get_m())
00151 {
00152
00153 if (movie)
00154 {
00155 movie_stream << "Current row: " << i << endl
00156 << "Find non-zero entry in the " << i
00157 << "-th row beginning from " << i << "-th entry\n";
00158 }
00159 j = i;
00160 while (j < m)
00161 {
00162 if (q(i, j) != 0) break;
00163 j++;
00164 }
00165 if (j >= m)
00166 {
00167
00168 q.del_row(i);
00169 e.ins_row(e.get_m(), f.take_row(i));
00170 if (movie)
00171 {
00172 movie_stream << "Zero row \n"
00173 << "Matrix is not of full rank \n";
00174 movie_print3();
00175 }
00176 continue;
00177 }
00178
00179 if (i != j)
00180 {
00181 q.swap_cols(i, j);
00182 a.swap_rows(i, j);
00183
00184 if (movie)
00185 {
00186 movie_stream << "Swap " << i << "-th and " << j << "-th columns:\n";
00187 movie_print3();
00188 }
00189 }
00190
00191 if (q(i, i) < 0)
00192 {
00193 q.mult_row(i, -1);
00194 f.mult_row(i, -1);
00195 }
00196
00197
00198 int_item b = q(i, i);
00199 index ii;
00200 for (ii = 0; ii < q.get_m(); ii++)
00201 if (ii != i)
00202 {
00203 int_item b_ii = q(ii, i);
00204 int_item alpha = gcd(b, b_ii);
00205 int_item b_i = b / alpha;
00206 b_ii = -b_ii / alpha;
00207 q.mult_row(ii, b_i);
00208 q.add_row(ii, i, b_ii);
00209 f.mult_row(ii, b_i);
00210 f.add_row(ii, i, b_ii);
00211
00212 alpha = gcd(gcd(q.row(ii)), gcd(f.row(ii)));
00213 q.div_row(ii, alpha);
00214 f.div_row(ii, alpha);
00215 }
00216 if (movie)
00217 {
00218 movie_stream << "Gaussian elimination in the " << i << "-th column: \n";
00219 movie_print3();
00220 }
00221 i++;
00222 }
00223 }
00224
00225
00226
00227
00228
00229
00230
00231
00232 template <class int_item>
00233 index motzkin_burger_algorithm<int_item>::fun(index j)
00234 {
00235
00236 index plus = 0;
00237 index minus = 0;
00238
00239
00240 for(index i = 0; i<current_ost; i++)
00241 {
00242 if(q(i, j)>0)
00243 plus++;
00244 if(q(i, j)<0)
00245 minus++;
00246 }
00247 return plus * minus;
00248 }
00249
00250
00251
00252
00253
00254
00255
00256 template <class int_item>
00257 index motzkin_burger_algorithm<int_item>::min_modification()
00258 {
00259 index min_column, jj, min_fun, cur_fun;
00260
00261
00262
00263
00264 min_column = current_column;
00265 min_fun = fun(min_column);
00266 for(jj = current_column + 1; jj < m; jj++)
00267 {
00268 cur_fun = fun(jj);
00269 if(cur_fun < min_fun)
00270 {
00271 min_column = jj;
00272 min_fun = cur_fun;
00273 }
00274 }
00275 return min_column;
00276 }
00277
00278
00279
00280
00281
00282
00283
00284 template <class int_item>
00285 index motzkin_burger_algorithm<int_item>::max_modification()
00286 {
00287 index max_column, jj, max_fun, cur_fun;
00288
00289
00290
00291
00292 max_column = current_column;
00293 max_fun = fun(max_column);
00294 for(jj = current_column + 1; jj < m; jj++)
00295 {
00296 cur_fun = fun(jj);
00297 if(cur_fun > max_fun)
00298 {
00299 max_column = jj;
00300 max_fun = cur_fun;
00301 }
00302 }
00303 return max_column;
00304 }
00305
00306
00307
00308
00309
00310
00311
00312
00313 template <class int_item>
00314 int motzkin_burger_algorithm<int_item>::reroof(index row)
00315 {
00316 for(index i = 0; i < common_zero.get_n(); i++)
00317 {
00318 index j = common_zero[i];
00319 if(q(row, j) != 0)
00320 return 0;
00321 }
00322 return 1;
00323 }
00324
00325
00326
00327
00328
00329
00330
00331 template <class int_item>
00332 int motzkin_burger_algorithm<int_item>::balance(index j_plus, index j_minus)
00333 {
00334
00335 for(index j = 0; j < current_column; j++)
00336 if(q(j_plus, j) == 0 && q(j_minus, j) == 0)
00337 common_zero.join(j);
00338
00339 if(common_zero.get_n()< r - 2)
00340
00341 return 0;
00342
00343
00344 for(index row = 0; row < current_ost; row++)
00345 if((row != j_plus) &&(row != j_minus))
00346 if(reroof(row))
00347 return 0;
00348
00349 return 1;
00350 }
00351
00352
00353
00354
00355
00356
00357
00358 template <class int_item>
00359 inline void motzkin_burger_algorithm<int_item>::movie_print3()
00360 {
00361 print3(movie_stream, transpose(a), f, q);
00362 movie_stream << endl;
00363 }
00364
00365
00366
00367
00368
00369
00370
00371 template <class int_item>
00372 void motzkin_burger_algorithm<int_item>::run()
00373
00374
00375
00376
00377
00378 {
00379 gauss();
00380
00381 if (movie)
00382 {
00383 movie_stream << "2. Motzkin algorithm\n"
00384 << " *****************\n";
00385 movie_print3();
00386 }
00387
00388 n = a.get_n();
00389 m = a.get_m();
00390 r = n - e.get_m();
00391 current_column = r;
00392 while(current_column < a.get_m())
00393 {
00394 vector<index> j_minus, j_plus;
00395
00396
00397
00398 current_ost = q.get_m();
00399
00400 index new_column;
00401 switch(modification)
00402 {
00403 case mm_no_modification:
00404 new_column = current_column;
00405 break;
00406 case mm_min_modification:
00407 new_column = min_modification();
00408 break;
00409 case mm_max_modification:
00410 new_column = max_modification();
00411 break;
00412 }
00413
00414 if (movie)
00415 {
00416 movie_stream << "Current column: " << current_column << endl;
00417 movie_stream << "New column: " << new_column << endl;
00418 }
00419
00420 if(current_column != new_column)
00421 {
00422 q.swap_cols(current_column, new_column);
00423 a.swap_rows(current_column, new_column);
00424
00425 if (movie)
00426 {
00427 movie_stream << "We swap columns: \n";
00428 movie_print3();
00429 }
00430 }
00431
00432
00433 for(index ii = 0; ii < current_ost; ii++)
00434 {
00435 if(q(ii, current_column) < 0)
00436 j_minus.join(ii);
00437 else if(q(ii, current_column) > 0)
00438 j_plus.join(ii);
00439 }
00440 if(j_minus.get_n() == current_ost)
00441 {
00442
00443 if (movie)
00444 movie_stream << "No solution exept zero\n";
00445 f.resize(0, n);
00446 e.resize(0, n);
00447 q.resize(0, m);
00448 return;
00449 }
00450 if(j_minus.get_n() == 0)
00451 {
00452
00453 if (movie)
00454 movie_stream << "This inequality is a corollary\n";
00455 a.del_row(current_column);
00456 q.del_col(current_column);
00457 m--;
00458 continue;
00459
00460 }
00461
00462 if (movie)
00463 movie_stream << "Balanced rows: ";
00464
00465 index jj_m;
00466 for(jj_m = 0; jj_m < j_minus.get_n(); jj_m++)
00467 {
00468 index j_m = j_minus[jj_m];
00469 for(index jj_p = 0; jj_p < j_plus.get_n(); jj_p++)
00470 {
00471 index j_p = j_plus[jj_p];
00472 common_zero.resize(0);
00473 if(balance(j_p, j_m))
00474 {
00475 if (movie)
00476 movie_stream << "(" << j_p << ", " << j_m << ") ";
00477 int_item b_minus = q(j_m, current_column);
00478 int_item b_plus = q(j_p, current_column);
00479 int_item delta = gcd(b_minus, b_plus);
00480 int_item alpha_plus = -b_minus / delta;
00481 int_item alpha_minus = b_plus / delta;
00482 f.add_rows(j_p, alpha_plus, j_m, alpha_minus);
00483 q.add_rows(j_p, alpha_plus, j_m, alpha_minus);
00484 delta = gcd(gcd(f.row(f.get_m() - 1)), gcd(q.row(q.get_m() - 1)));
00485 f.div_row(f.get_m()-1, delta);
00486 q.div_row(q.get_m()-1, delta);
00487 }
00488 }
00489 }
00490
00491 if (movie)
00492 {
00493 movie_stream << endl;
00494 movie_print3();
00495 movie_stream << "Delete negative rows:\n";
00496 }
00497
00498 f.del_rows(j_minus);
00499 q.del_rows(j_minus);
00500
00501 if (movie) movie_print3();
00502
00503 current_column++;
00504 }
00505 }
00506
00507 #endif