API tools faq
paste
Advertisement
hugol

Untitled

hugol
Apr 27th, 2015
209
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 11.73 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <vector>
  3. #include <assert.h>
  4. #include <math.h>
  5. #include <map>
  6. #include <time.h>
  7.  
  8. using namespace std;
  9.  
  10. #define d 0.85
  11. #define MAX_ITER 1000
  12. #define EPSILON 0.0000001
  13.  
  14. #define ZERO 1e-12 // stała przybliżenia zera
  15.  
  16.  
  17. void setXY(map<int, map<int, double>> &macierz_sasiedztwa, int x, int y, double val)
  18. {
  19.     if (fabs(val) >= ZERO)
  20.         macierz_sasiedztwa[x][y]= val;
  21. }
  22.  
  23. double getXY(map<int, map<int, double>> &macierz_sasiedztwa, int x, int y)
  24. {
  25.     if ((macierz_sasiedztwa.find(x) != macierz_sasiedztwa.end())
  26.         && (macierz_sasiedztwa[x].find(y) != macierz_sasiedztwa[x].end()) )
  27.     {
  28.         return macierz_sasiedztwa[x][y];
  29.     }
  30.     else
  31.         return 0;
  32.  
  33. }
  34.  
  35. double getXYRazyMacierzDiag(vector <double> &ilosc, map<int, map<int, double>> &macierz_sasiedztwa, int x, int y)
  36. {
  37.         if (y>=ilosc.size())
  38.                 return getXY(macierz_sasiedztwa, x, y);
  39.  
  40.  
  41.         double ret = 0;
  42.         if ((macierz_sasiedztwa.find(x) != macierz_sasiedztwa.end())
  43.                 && (macierz_sasiedztwa[x].find(y) != macierz_sasiedztwa[x].end()) )
  44.         {
  45.                 // jezeli jest polaczenie
  46.                 ret = ret - macierz_sasiedztwa[x][y];
  47.                 ret = ret * (double)d/(double)ilosc[y];
  48.                 ret = ret + (x==y ? 1.0 : 0.0);
  49.                 // return -(double)macierz_sasiedztwa[x][y] * (double)d/(double)ilosc[y] + x==y;
  50.         }
  51.         else
  52.         {
  53.                 if (!(ilosc[x] > 0))
  54.                 {
  55.                                 //return double(0.0-(1.0 * d/(double)ilosc.size()) + x==y);
  56.                                 ret = ret - (1.0 * d/(double)ilosc.size());
  57.                                 ret = ret + (double)x==y;
  58.                 }
  59.                 else
  60.                         //return double(0.0 + x==y);
  61.                         ret = (double)x==y;
  62.         }
  63.  
  64.         return ret;
  65.  
  66. }
  67.  
  68. void setXYRazyMacierzDiag(vector <double> &ilosc, map<int, map<int, double>> &macierz_sasiedztwa, int x, int y, double val)
  69. {
  70.         if (fabs(val) < ZERO)
  71.                 return;
  72.  
  73.  
  74.         if (y>=ilosc.size())
  75.         {
  76.                 setXY(macierz_sasiedztwa, x, y, val);
  77.                 return;
  78.         }
  79.  
  80.         double set;
  81.  
  82.         // jezeli jest polaczenie
  83.         //ret = ret - macierz_sasiedztwa[x][y];
  84.         //ret = ret * (double)d/(double)ilosc[y];
  85.         //ret = ret + (x==y ? 1.0 : 0.0);
  86.  
  87.         set = val;
  88.         set = set - (x==y ? 1.0 : 0.0);
  89.         set = set / ((double)d/(double)ilosc[y]);
  90.         macierz_sasiedztwa[x][y] = -set;
  91.  
  92. }
  93. // sprawdzanie
  94. bool sprawdz_zbieznosc(vector <vector<double>> in)
  95. {
  96.     assert( in.size() == in[0].size());
  97.  
  98.     // ||C|| W
  99.     double max = 0;
  100.     for(int i=0; i<in.size(); i++)
  101.     {
  102.         double sum = 0;
  103.         for(int j=0; j<in.size(); j++)
  104.         {
  105.             sum += abs( in[i][j] );
  106.         }
  107.         if (sum > max)
  108.             max = sum;
  109.     }
  110.     if (max < 1)
  111.         return true;
  112.  
  113.     // ||C|| K
  114.     max = 0;
  115.     for(int i=0; i<in.size(); i++)
  116.     {
  117.         double sum = 0;
  118.         for(int j=0; j<in.size(); j++)
  119.         {
  120.             sum += abs( in[j][i] );
  121.         }
  122.         if (sum > max)
  123.             max = sum;
  124.     }
  125.     if (max < 1)
  126.         return true;
  127.  
  128.     // ||C|| F
  129.     double sum_sum = 0;
  130.     for(int i=0; i<in.size(); i++)
  131.     {
  132.         double sum = 0;
  133.         for(int j=0; j<in.size(); j++)
  134.         {
  135.             sum += abs( in[j][i]*in[j][i] );
  136.         }
  137.         sum_sum + sum;
  138.     }
  139.     if (sqrt(sum_sum) < 1)
  140.         return true;
  141.  
  142.     return false;
  143. }
  144.  
  145. void odwroc_macierz(vector<vector<double>> &a)
  146. {
  147.     int i,j,k,n=a.size();
  148.     double d_tmp;        
  149.     // rozszerzenie o macierz diagonalna
  150.     for(i=0;i<n;i++)
  151.     {
  152.         for(j=0;j<n;j++)
  153.         {
  154.             if(j==i)
  155.                 a[i].push_back( 1 );
  156.             else
  157.                 a[i].push_back( 0 );
  158.         }
  159.     }
  160.  
  161.     /************** partial pivoting **************/
  162.     for(i=n-1;i>0;i--)
  163.     {
  164.         if(a[i-1][0]<a[i][0])
  165.         {
  166.             for(j=0;j<n*2;j++)
  167.             {
  168.                 d_tmp=a[i][j];
  169.                 a[i][j]=a[i-1][j];
  170.                 a[i-1][j]=d_tmp;
  171.             }
  172.         }
  173.     }
  174.  
  175.     /********** reducing to diagonal  matrix ***********/
  176.  
  177.     for(i=0;i<n;i++)
  178.     {
  179.         for(j=0;j<n;j++)
  180.             if(j!=i)
  181.             {
  182.                 d_tmp=a[j][i]/a[i][i];
  183.                 for(k=0;k<n*2;k++)
  184.                     a[j][k]-=a[i][k]*d_tmp;
  185.             }
  186.     }
  187.     /************** reducing to unit matrix *************/
  188.     for(i=0;i<n;i++)
  189.     {
  190.         d_tmp=a[i][i];
  191.         for(j=0;j<n*2;j++)
  192.             a[i][j]=a[i][j]/d_tmp;
  193.     }
  194.  
  195.     for(i=0;i<n;i++)
  196.     {
  197.         for(j=0;j<n;j++)
  198.             a[i].erase( a[i].begin() );
  199.     }
  200.  
  201. }
  202.  
  203. void odejmij_macierz(vector<vector<double>> &A, vector<vector<double>> &B)
  204. {
  205.     assert(A.size() == B.size());
  206.     assert(A[0].size() == B[0].size());
  207.  
  208.     for(int i=0; i<A.size(); i++)
  209.     {
  210.         for(int j=0; j<A[0].size(); j++)
  211.         {
  212.             A[i][j] -= B[i][j];
  213.         }
  214.     }
  215. }
  216.  
  217. void pomnoz_macierz(vector<vector<double>> &A, vector<vector<double>> &B)
  218. {
  219.     for(int i = 0; i < A.size(); i++)
  220.     {
  221.         for(int j = 0; j < B[0].size(); j++)
  222.         {
  223.             double s = 0;
  224.             for(int k = 0; k < A[0].size(); k++)
  225.                 s += A[i][k] * B[k][j];
  226.             A[i][j] = s;
  227.         }
  228.     }
  229. }
  230.  
  231. int Gauss_Seidel(map<int, map<int, double>> &A, vector<double> &B,  vector<double> &X, vector <double> &ilosc)
  232. {
  233.     // macierz kwadratowa musi byc
  234.     //assert(A.size() == A[0].size());
  235.     //assert(A.size() == B.size());
  236.  
  237.     // Calculate D^-1
  238.     for (map<int, map<int, double>>::iterator it=A.begin(); it!=A.end(); ++it)
  239.     {
  240.         for(map<int, double>::iterator it2=it->second.begin(); it2!=it->second.end(); ++it2)
  241.         {
  242.             if (it->first == it2->first)
  243.                 setXYRazyMacierzDiag(ilosc, A, it->first, it2->first, 1.0 / getXYRazyMacierzDiag(ilosc, A, it->first, it2->first));
  244.         }
  245.     }
  246.  
  247.     // Calculate D^-1 * b
  248.     for (map<int, map<int, double>>::iterator it=A.begin(); it!=A.end(); ++it)
  249.     {
  250.         int i = it->first;
  251.         B[i] *= getXYRazyMacierzDiag(ilosc, A, i, i);
  252.     }
  253.  
  254.     //Calculate D^-1 * L
  255.     //Calculate D^-1 * U
  256.     for (map<int, map<int, double>>::iterator it=A.begin(); it!=A.end(); ++it)
  257.     {
  258.         for(map<int, double>::iterator it2=it->second.begin(); it2!=it->second.end(); ++it2)
  259.         {
  260.             int i = it->first;
  261.             int j = it2->first;
  262.             if (j == i)
  263.                 continue;
  264.             //A[i][j] *= A[i][i];
  265.             setXYRazyMacierzDiag(ilosc, A, i, j, getXYRazyMacierzDiag(ilosc, A, i, j) *  getXYRazyMacierzDiag(ilosc, A, i, i) );
  266.         }
  267.     }
  268.  
  269.  
  270.     // iteracje
  271.     int k;
  272.     for (k=0; k<MAX_ITER; k++)
  273.     {
  274.         double maxXiDiff = 0, maxXi = 0;
  275.         double X2;
  276.         //for (int i=0; i<A.size(); i++)
  277.         for (map<int, map<int, double>>::iterator it=A.begin(); it!=A.end(); ++it)
  278.         {
  279.             int i = it->first;
  280.             X2 = B[i];  
  281.             // x = D^-1*b -
  282.             for(map<int, double>::iterator it2=it->second.begin(); it2!=it->second.end(); ++it2)
  283.             {
  284.                 int j = it2->first;
  285.                 if (j == i)
  286.                     continue;
  287.                 //X2 -= A[i][j]*X[j];    // - D^-1*L * x
  288.                 X2 -= getXYRazyMacierzDiag(ilosc, A, i, j)*X[j];    // - D^-1*L * x
  289.             }
  290.             if ( abs( X2 - X[i] ) > maxXiDiff)
  291.                 maxXiDiff=abs( X2 - X[i] );
  292.             if ( abs( X2 ) > maxXi )
  293.                 maxXi = abs( X2 );
  294.             X[i] = X2;
  295.         }
  296.         double eps = maxXiDiff / maxXi;
  297.         if ( !( EPSILON < eps) )
  298.             break;
  299.     }
  300.  
  301.     return k;
  302. }
  303.  
  304. int CarlGustavJakobJacobi(vector<vector<double>> &A, vector<double> &B,  vector<double> &X)
  305. {
  306.     // macierz kwadratowa musi byc
  307.     assert(A.size() == A[0].size());
  308.     assert(A.size() == B.size());
  309.  
  310.     // wyznaczenie macierzy N = D^-1
  311.     vector<double> N;
  312.     for (int i=0; i<A.size(); i++)
  313.     {
  314.         N.push_back( 1.0/A[i][i] );
  315.     }
  316.  
  317.     // kryt zbieznosci wedlug wykladu
  318.  
  319.     // wyznaczenie macierzy M = -N (L + U)
  320.     vector<vector<double>> M;
  321.     for (int i=0; i<A.size(); i++)
  322.     {
  323.         M.push_back( vector<double>() );
  324.         for (int j=0; j<A.size(); j++)
  325.         {
  326.             if (i == j)
  327.             {
  328.                 M[i].push_back( 0 );
  329.             }
  330.             else
  331.             {
  332.                 M[i].push_back ( -(A[i][j] * N[i]) );
  333.             }
  334.         }
  335.     }
  336.     assert( sprawdz_zbieznosc(M) );
  337.  
  338.     vector<double> X2 = X;
  339.     int k; // ilosc iteracji
  340.     for (k=0; k<MAX_ITER; k++)
  341.     {
  342.         double maxXiDiff = 0, maxXi = 0;
  343.         for (int i=0; i<A.size(); i++)
  344.         {
  345.             X2[i] = N[i]*B[i];
  346.             for (int j=0; j<A.size(); j++)
  347.             {
  348.                 X2[i] += M[i][j]*X[j];
  349.             }
  350.             if ( abs( X2[i] - X[i] ) > maxXiDiff)
  351.                 maxXiDiff = abs( X2[i] - X[i] );
  352.             if ( abs( X2[i] ) > maxXi )
  353.                 maxXi = abs( X2[i] );
  354.         }
  355.         X = X2;
  356.         double eps = maxXiDiff / maxXi;
  357.         if ( !( EPSILON < eps) )
  358.             break;
  359.     }
  360.  
  361.     return k;
  362. }
  363.  
  364. void generuj_wektor_poczatkowy( vector<double> &X, int size)
  365. {
  366.     for (int i=0; i<size; i++)
  367.     {
  368.         X.push_back( 0.0 );
  369.     }
  370. }
  371.  
  372. void PR_dodaj_kolumne(int n, vector<double> &B)
  373. {
  374.     for(int i=0; i<n; i++)
  375.     {
  376.         B.push_back( (1.0-d) / double(n) );
  377.     }
  378. }
  379.  
  380. void PR_pomnoz_razy_macierz_diagonalna_i_dodaj_krawedzie_do_wszystkich_pustych
  381.     (vector<double> &ilosc, map<int, map<int, double>> &sasiedztwo)
  382. {
  383.  
  384.     for(int i=0; i<ilosc.size(); i++)
  385.     {
  386.         // dodanie krawedzi do wszystkich
  387.         if (!(ilosc[i] > 0))
  388.         {
  389.             ilosc[i] = ilosc.size()-1;
  390.             for (int j=0; j<ilosc.size(); j++)
  391.             {
  392.                 if (i==j)
  393.                     continue;
  394.                 sasiedztwo[j][i] = 1;
  395.             }
  396.         }
  397.  
  398.  
  399.         //for(int j=0; j<sasiedztwo[i].size(); j++)
  400.         for(map<int, map<int, double>>::iterator it=sasiedztwo.begin(); it!=sasiedztwo.end(); ++it)
  401.         {
  402.             if (it->second.find(i) != it->second.end() )
  403.                 it->second[i] *= d/ilosc[i];
  404.         }
  405.     }
  406. }
  407.  
  408. void PR_I_odj_macierz(map<int, map<int, double>> &sasiedztwo, int n)
  409. {
  410.     for (map<int, map<int, double>>::iterator it=sasiedztwo.begin(); it!=sasiedztwo.end(); ++it)
  411.     {
  412.         for(map<int, double>::iterator it2=it->second.begin(); it2!=it->second.end(); ++it2)
  413.         {
  414.             double val = 0 - getXY(sasiedztwo, it->first, it2->first);
  415.             val += it->first==it2->first;
  416.             setXY(sasiedztwo, it->first, it2->first, val);
  417.         }
  418.     }
  419.     for(int i=0; i<n; i++)
  420.     {
  421.         if (getXY(sasiedztwo, i, i) == 0)
  422.             setXY(sasiedztwo, i, i, 1);
  423.     }
  424. }
  425.  
  426. bool PR_wczytaj(map<int, int> &backNodes, vector <double> &ilosc, map<int, map<int, double> > &macierz_sasiedztwa)
  427. {
  428.     int currId = 0;
  429.  
  430.     int src, dst;
  431.     int count = 0;
  432.  
  433.     map<int, int> Nodes;
  434.  
  435.     while (!cin.eof())
  436.     {
  437.         cin >> src >> dst;
  438.         if (cin.eof() || src == -1)
  439.             break;
  440.  
  441.         // mapy id
  442.         if (Nodes.find(src) == Nodes.end())
  443.         {
  444.             Nodes[src] = currId;
  445.             backNodes[currId] = src;
  446.             currId++;
  447.         }
  448.         if (Nodes.find(dst) == Nodes.end())
  449.         {
  450.             Nodes[dst] = currId;
  451.             backNodes[currId] = dst;
  452.             currId++;
  453.         }
  454.  
  455.         /* powiekszanie wektorow*/
  456.         if (Nodes[src] > count)
  457.             count = Nodes[src];
  458.         if (Nodes[dst] > count)
  459.             count = Nodes[dst];
  460.         while(ilosc.size() <= count)
  461.             ilosc.push_back( 0 );
  462.  
  463.  
  464.         // ustawienie macierzy
  465.         ilosc[Nodes[src]]++;
  466.         macierz_sasiedztwa[Nodes[dst]][Nodes[src]] = 1;
  467.  
  468.     }
  469.  
  470.     return true;
  471. }
  472.  
  473. #define PAGE_RANK
  474.  
  475. int main()
  476. {
  477.     map<int, map<int, double> > A;
  478.     vector<double> B;
  479.     vector<double> X;
  480.  
  481.  
  482. #ifndef PAGE_RANK
  483.     int n;
  484.     cin >> n;
  485.     // wczytywanie macierzy
  486.     for(int i=0; i<n; i++)
  487.     {
  488.         double tmp;
  489.         for(int j=0; j<n; j++)
  490.         {
  491.             cin >> tmp;
  492.             A[i][j] = tmp;
  493.         }
  494.         cin >> tmp;
  495.         B.push_back(tmp);
  496.     }
  497.     generuj_wektor_poczatkowy(X, n);
  498. #endif
  499.  
  500.  
  501. #ifdef PAGE_RANK
  502.     vector<double> ilosc;
  503.     map<int, int> Nodes;
  504.     PR_wczytaj(Nodes, ilosc, A);   
  505.  
  506.     /*PR_pomnoz_razy_macierz_diagonalna_i_dodaj_krawedzie_do_wszystkich_pustych(ilosc, A);
  507.     PR_I_odj_macierz(A, ilosc.size() );*/
  508.  
  509.     PR_dodaj_kolumne(ilosc.size(), B);
  510.     generuj_wektor_poczatkowy(X, ilosc.size() );
  511. #endif
  512.  
  513.     time_t now; time(&now);
  514.     int iter = Gauss_Seidel(A, B, X, ilosc);
  515.     time_t seconds = time(NULL) - now;
  516.     cout << "Gauss_Seidel: liczba iteracji: " << iter << endl;
  517.     cout << "Czas [s]: " << seconds << endl;
  518.     double sum = 0;
  519.     for(int i=0; i<X.size(); i++)
  520.     {
  521. //#ifndef PAGE_RANK
  522.         cout << "x" << i << "= " << X[i] << endl;
  523. //#endif
  524.         sum += X[i];
  525.     }
  526.     cout << "Suma: " << sum << endl;
  527.  
  528.     X.clear();
  529.     generuj_wektor_poczatkowy(X, ilosc.size() );
  530.  
  531.     time(&now);
  532.     /*iter = CarlGustavJakobJacobi(A, B, X);*/
  533.     seconds = time(NULL) - now;
  534.     cout << "CarlGustavJakobJacobi: liczba iteracji: " << iter << endl;
  535.     cout << "Czas [s]: " << seconds << endl;
  536.     sum = 0;
  537.     for(int i=0; i<X.size(); i++)
  538.     {
  539. #ifndef PAGE_RANK
  540.         cout << "x" << i << "= " << X[i] << endl;
  541. #endif
  542.         sum += X[i];
  543.     }
  544.     cout << "Suma: " << sum << endl;
  545.  
  546.     return 0;
  547. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!