Zobrazení zdroje stránky GNU Gama, vyrovnání měření zprostředkujících s diagonální maticí vah

Součástí projektu [[GNU Gama]] je C++ knihovna tříd a funkcí. Třída <code>Adj</code> implementuje úlohu ''vyrovnání měření zprostředkujících''

<center>
<math>
\mathbf{A}\mathbf{x} = \mathbf{l}, \qquad \mathbf{P}
</math>
</center>

kde '''A''' je první matice plánu, '''x''' vektor určovaných parametrů, '''l''' vektor absolutních členů a '''P''' matice vah.

Pro řídkou matici '''A''' a diagonální váhovou matici '''P'''
může zapsat vstupní data v následujícím formátu

   ''počet_neznámých   počet_měření''

následuje popis jednotlivých řádků řídké soustavy

   ''N  index_1  index_2 ... index_N''
   ''váha  pravá_strana  koef_1  koef_2  ... koef_N''

Například soustavu

<center>
<math>
\mathbf{A} =  \begin{pmatrix}  
   1.0 &   0  &   0  &   0  \\
   0.9 &   0  &  0.1 &   0  \\ 
   2.0 &   0  &  3.0 &  0.1 \\ 
   0.5 &  0.6 &  0.2 &  0.3 \\
  -0.1 &   0  &  0.1 &  3.0 \\
   0.6 & -0.1 & -0.9 &  0.4 
\end{pmatrix},\qquad
\mathbf{l} = \begin{pmatrix} 
    1.03 \\
    1.18 \\
   11.38 \\
    3.51 \\
   12.21 \\
   -0.66 
\end{pmatrix},\qquad
\mathbf{P} = \mathrm{diag}\begin{pmatrix} 
   1  \\
  1.1 \\
  1.2 \\
  1.3 \\
  1.4 \\
  1.5 
\end{pmatrix}
</math>
</center>

zapíšeme jako 

<pre>
   4  6

   1          1
   1   1.03   1
   2          1    3
  1.1  1.18  0.9  0.1
   3          1    3    4
  1.2 11.38  2.0  3.0  0.1   
   4          1    2    3    4
  1.3  3.51  0.5  0.6  0.2  0.3
   3          1    3    4
  1.4 12.21 -0.1  0.1  3.0
   4          1    2    3    4
  1.5 -0.66  0.6 -0.1 -0.9  0.4
</pre>

počet mezer není významný, zde jsou vícenásobné mezery použity pouze pro formátování pro lepší čitelnost.  Tento vstupní formát čte metoda <code>read_gama_local_old_format(cin)</code> třídy
<code>AdjInputData</code>.

Požití obou tříd ukazuje následující demo program
<pre>
/*******************************************************************

Vstupni data:
------------

Zapisuji se pouze nenulove koeficienty rovnic oprav. Prvni zaznam
obsahuhue pouze dve cela cisla

   pocet_neznamych  pocet_mereni

dale pak nasleduji dva radky pro kazdou rovnici

   N index_1 index_2 ... index_N
   vaha prava_strana koef_1 koef_2 ... koef_N

kde N je pocet nenulovych prvku v danem radku


Volani programu:
---------------

   kombinace-01 [envelope|gso|svd|cholesky] < vstup.txt > vystup.xml

Implicitni algoritmus je reseni ridke matice v pametovem modelu
envelope.

********************************************************************/

#include <iostream>
#include <sstream>
#include <cmath>
#include <vector>
#include <gnu_gama/adj/adj.h>

using std::cin;
using std::cout;
using std::endl;
using std::cerr;
using std::sqrt;
using std::vector;
using std::string;


int main(int argc, char* argv[])
{
  using namespace GNU_gama;

  // objekt typu AdjInputData musi byt vytvoren dynamicky, protoze
  // pointer prejima objekt Adj, ktery pro vtupni data vola destruktor

  AdjInputData* data = AdjInputData();
  data->read_gama_local_old_format(cin);

  Adj adj;
  adj.set(data);

  string alg = "envelope";
  if (argc == 2)
    {
      if      (alg == string("envelope")) adj.set_algorithm(Adj::envelope);
      else if (alg == string("gso")     ) adj.set_algorithm(Adj::gso);
      else if (alg == string("svd")     ) adj.set_algorithm(Adj::svd);
      else if (alg == string("cholesky")) adj.set_algorithm(Adj::cholesky);
      else
        {
          cerr <<
            "volani : prog [envelope|gso|svd|cholesky] "
            "< vstup > vystup\n\n";
        }
    }

  Vec<> x = adj.x();
  Vec<> r = adj.r();

  cout << "# algoritmus = " << alg << endl;
  cout << "# nezname    = " << x.dim() << endl;
  cout << "# opravy     = " << r.dim() << endl;
  cout << "# defekt     = " << adj.defect() << endl;
  cout << "# [pvv]      = " << adj.rtr() << endl;
  cout << endl;

  double m02 = r.dim() > x.dim() ? adj.rtr()/(r.dim() - x.dim()) : 0.0;

  for (int i=1; i<=x.dim(); i++)
    {
      cout << i << "\t" << x(i) << "\t" << sqrt(m02*adj.q_xx(i,i)) << endl;
    }

  // cout << endl;
  //
  // for (int i=1; i<=r.dim(); i++)
  //   {
  //     cout << i << "\t" << r(i) << "\t" << sqrt(m02*adj.q_bb(i,i)) << endl;
  //   }
}
</pre>

Výstup z programu pro uvedený příklad je

<pre>
# algoritmus = envelope
# nezname    = 4
# opravy     = 6
# defekt     = 0
# [pvv]      = 0.00158437

1       1.01759 0.0142184
2       2.00305 0.0410978
3       2.97991 0.0111895
4       4.00515 0.00783747
</pre>

[[Kategorie:GNU]]