Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members   Examples  

motzkin.h

Go to the documentation of this file.
00001 //N.Yu.Zolotykh 1999
00002 //University of Nizhni Novgorod, Russia
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  I m p l e m e n t a t i o n
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 /*          gauss          */
00114 /*                         */
00115 /***************************/
00116 
00117 /*Элементарными преобразованиями над строками матрицы a
00118   и перестановкой ее столбцов приводит ее к диагональному виду q:
00119     f * transpose(a) = q.*/
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   // error handling?
00131 
00132   /* f = identity */
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   /* main loop */
00149   i = 0;
00150   while (i < q.get_m())
00151   {
00152     /* find non-zero entry in the i-th row beginning from i-th entry */
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       /* it's a zero row */
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; // main loop
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     /* В i-ом столбце образуем все нули */
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++; //next row
00222   }
00223 }
00224 
00225 
00226 /***************************/
00227 /*                         */
00228 /*           fun           */
00229 /*                         */
00230 /***************************/
00231 
00232 template <class int_item>
00233   index motzkin_burger_algorithm<int_item>::fun(index j)
00234 {
00235   //must be indexindex ! !!!!!
00236   index plus = 0;
00237   index minus = 0;
00238   //the number of positive and negative entries in the j - th column
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 /*    min_modification     */
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   /* the column with minimal product of the number of positive
00262      and the number of negative entries */
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 /*    max_modification     */
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   /* the column with maximal product of the number of positive
00290      and the number of negative entries */
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 /*         reroof          */
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 /*         balance         */
00328 /*                         */
00329 /***************************/
00330 
00331 template <class int_item>
00332   int motzkin_burger_algorithm<int_item>::balance(index j_plus, index j_minus)
00333 {
00334   //fist we construct the set CommonZero
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 /*       movie_print3      */
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 /*           run           */
00368 /*                         */
00369 /***************************/
00370 
00371 template <class int_item>
00372   void motzkin_burger_algorithm<int_item>::run()
00373      /*
00374        find extreme rays f and basis e of a cone ax>= 0
00375        returns also an incidence matrix q = f * transpose(a)
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     //for each inequality from a, beginning from(r + 1) th
00397       /* show_algorithm_state("current_column = ", current_column); */
00398     current_ost = q.get_m();
00399     //taking a new inequality depends upon modification
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     //now we are constructing the sets j_plus, j_minus
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       //there is no solution except 0
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       //this is a Н еравенство-следствие
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       //new inequality
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       } //Iterate J_plus
00489     } //Iterate J_minus
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

Generated at Tue Jan 22 20:37:04 2002 for Arageli by doxygen1.2.9.1 written by Dimitri van Heesch, © 1997-2001