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 Adj implementuje úlohu vyrovnání měření zprostředkujících

$\mathbf{A}\mathbf{x} = \mathbf{l}, \qquad \mathbf{P}$

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

$\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}$

zapíšeme jako

   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

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 read_gama_local_old_format(cin) třídy AdjInputData.

Požití obou tříd ukazuje následující demo program

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

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;
  //   }
}

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

# 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