Aproksymacja średniokwadratowa punktowa wielomianowa

/*
Dawid Mocek

Działa pod VS 2013 lub kompilatorem z C++11
Ten kod w 100% działa pod windows
Program wygeneruje od groma plików - lepiej je poprzenosic w osobne katalogi. Nazwy plików są schematyczne

Aproksymacja średniokwadratowa punktowa wielomianowa funkcji:
y = sinx * cos3x

Ogólna treść zadania:
Przeprowadzić aproksymację średniokwadratową punktowąw ielomianową
funkcji y = sinx * cos3x x w przedziale <-2,2> dla 50 węzłów.
Zwiększyć stopień wielomianu aproksymacyjnego od 1 do 20 i szukać najmniejszego błędu aproksymacji.
Zrobić również doświadczenia na innych wartościach.

Co trzeba zmienic:
string katalog_wynikowy

Co można zmienic:

DELIMITER
DIR_DELIMITER

string katalog_wynikowy
int liczbawezlow
vector<int> stopienie
vector<vector<double>> granice

Pliki CSV można otworzyć bezpośrednio w Excelu
*/

#include <iostream>
#include <cmath>
#include <fstream>
#include <string>
#include <vector>
#include <sstream>
#include <algorithm>

// Excel potrzebuje dla naszego jezyka(locale) takiego rozdzielacza dla pików csv:
#define DELIMITER ';'

// Windows - pod Linuksem stosować "/"
#define DIR_DELIMITER "\\"

using namespace std;

// Dla czytelniejszego kodu
typedef vector<double> Wektor;
typedef vector<vector<double>> Macierz;

// Struktura przechowuje dane-logi po każdej iteracji
// Zastosowanie: łatwe znalezienie najmniejszego błędu apro. w plikach
// raportów
struct DaneLog
{
    double blad_aproksymacji;
    double granica_prawa;
    double granica_lewa;
    int stopien;
    string plik_raportu;
    string plik_excela;

    DaneLog();
    DaneLog(double blad_aproksymacji, double granica_prawa, double granica_lewa, int stopien, string plik_raportu, string plik_excela);
    bool operator< (const DaneLog& org) const;
    bool operator> (const DaneLog& org) const;
    bool operator== (const DaneLog& org) const;
    bool operator!= (const DaneLog& org) const;
    friend ostream& operator<< (ostream& out, const DaneLog& org);
};

DaneLog::DaneLog()
: blad_aproksymacji(0.0), granica_prawa(0.0), granica_lewa(0.0), stopien(0), plik_excela(""), plik_raportu("")
{

};

DaneLog::DaneLog(double blad_aproksymacji, double granica_prawa, double granica_lewa, int stopien, string plik_raportu, string plik_excela)
: blad_aproksymacji(blad_aproksymacji), granica_prawa(granica_prawa), granica_lewa(granica_lewa), stopien(stopien), plik_raportu(plik_raportu), plik_excela(plik_excela)
{

};

bool DaneLog::operator< (const DaneLog& org) const
{
    return blad_aproksymacji < org.blad_aproksymacji;
};
bool DaneLog::operator> (const DaneLog& org) const
{
    return !DaneLog::operator<(org);
};
bool DaneLog::operator== (const DaneLog& org) const
{
    return blad_aproksymacji == org.blad_aproksymacji;
};
bool DaneLog::operator!= (const DaneLog& org) const
{
    return !DaneLog::operator==(org);
};
ostream& operator<< (ostream& out, const DaneLog& org)
{
    return out << "Blad (st. = " << org.stopien << " ) = " << org.blad_aproksymacji << " w raporcie: \"" << org.plik_raportu << "\" i Excelu: \"" << org.plik_excela << "\", dla g.l. " << org.granica_lewa
        << " i g.p. " << org.granica_prawa;
}

///////////////////////////////////////////////////////////////// Glowne funkcje /////////////////////////////////////////////////////////////////

// Liczy G
Macierz * _G(const Wektor *wezly, const int stopien, const int liczbawezlow)
{
    Macierz *G = new Macierz();
    G->resize(stopien + 1);
    for (int z = 0; z <= stopien; z++)
        (*G)[z].resize(stopien + 1);

    for (int j = 0; j <= stopien; ++j)
    for (int k = 0; k <= stopien; ++k)
    {
        (*G)[k][j] = 0;

        for (int i = 0; i < liczbawezlow; ++i)
            (*G)[k][j] += pow((*wezly)[i], k + j);

    }

    return G;
}

// Liczy R
Wektor * _R(const Wektor *wezly, const Wektor * fx, const int stopien, const int liczbawezlow)
{
    Wektor *R = new Wektor();
    R->resize(stopien + 1);
    for (int i = 0; i <= stopien; ++i)
    {
        (*R)[i] = 0;
        for (int j = 0; j < liczbawezlow; ++j)
            (*R)[i] += (*fx)[j] * pow((*wezly)[j], i);
    }
    return R;
}

// Przygotowuje wezly oraz wartość f(x)
void _Przygotuj(Wektor *fx, Wektor *wezly, const int liczbawezlow, const double granica_lewa, const double granica_prawa)
{
    fx->resize(liczbawezlow);
    wezly->resize(liczbawezlow);

    double tmp = (granica_prawa - granica_lewa) / (liczbawezlow - 1);
    for (int i = 0; i < liczbawezlow; ++i)
    {
        double x = granica_lewa + (i * tmp);
        // y = sinx * cos3x
        (*fx)[i] = sin(x)*(cos(3 * (x)));
        (*wezly)[i] = x;
    }
}

// Macierz trójątna - zwraca false jesli dojdzie do dzielenia przez zero
bool _MacierzTrj(Wektor *R, Macierz *G, const int stopien)
{

    for (int k = 0; k <= stopien; ++k)
    {
        for (int i = k + 1; i <= stopien; ++i)
        {
            (*R)[i] -= ((*G)[i][k] * (*R)[k]) / (*G)[k][k];
            for (int j = k + 1; j <= stopien; ++j)
            {
                if ((*G)[k][k] == 0) return false;
                (*G)[i][j] -= ((*G)[i][k] * (*G)[k][j]) / (*G)[k][k];
            }
            (*G)[i][k] = 0;
        }
    }
    return true;
}

// Wyznaczamy niewiadomoe
Wektor * _A(const Wektor *R, const Macierz *G, const int stopien)
{
    double suma;
    Wektor *A = new Wektor();
    A->resize(stopien + 1);

    (*A)[stopien] = (*R)[stopien] / (*G)[stopien][stopien];
    for (int i = stopien - 1; i >= 0; --i)
    {
        suma = 0.0;
        for (int j = i + 1; j <= stopien; ++j)
            suma += (*G)[i][j] * (*A)[j];

        (*A)[i] = ((*R)[i] - suma) / (*G)[i][i];
    }
    return A;
}

// Wyznacza wartości funkcji aproksymujacej
Wektor * _Qx(const Wektor *A, const Wektor *_wezly, const int stopien, const int liczbawezlow)
{
    Wektor *Qx = new Wektor();
    Qx->resize(liczbawezlow);
    for (int i = 0; i < liczbawezlow; ++i)
    for (int j = 0; j <= stopien; ++j)
        (*Qx)[i] += (*A)[j] * pow((*_wezly)[i], j);

    return Qx;
}

// Oblicza błąd aproksymacji
double _Blad(const Wektor *Qx, const Wektor *Fx, const int liczbawezlow)
{
    double tmp = 0;
    double blad_aproksymacji = 0;
    for (int i = 0; i < liczbawezlow; ++i)
        tmp += pow((*Qx)[i] - (*Fx)[i], 2);

    return (sqrt(tmp / liczbawezlow));

}
//
////////////////////////////////////// Koniec głównych funkcji //////////////////////////////

// Raport do pliku
string  _Raport(const string dir, const Wektor *wezly, const Wektor *fx, const Wektor *qx, const Wektor *A, const int liczbawezlow, const int stopien, const double granica_lewa, const double granica_prawa, const double blad)
{
    ostringstream  ss_file;
    ss_file << "Raport_Stopien_" << stopien << "-liczbaWezlow_" << liczbawezlow << "-granicaL_" << granica_lewa << "-granicaP_" << granica_prawa << ".txt";
    string plik = ss_file.str();
    fstream file(dir + (string)DIR_DELIMITER + plik, ios::out | ios::trunc);
    file.precision(20);
    file.imbue(locale(""));
    file << "+-------------------------------------------------------------------------------------------------------+" << endl;
    file << "| Przedzial: [" << granica_lewa << ", " << granica_prawa << "]" << endl;
    file << "| Stopien: " << stopien << endl;
    file << "| Liczba wezlow: " << liczbawezlow << endl;
//  file.setf(ios::scientific);
    file << "| Blad aproksymacji: " << blad << endl;

    file.setf(ios::right, ios::adjustfield);
    file.setf(ios::fixed, ios::floatfield);
    file << "| Tabela wartości funkcji stablicowanej i aproksymujacej w wezlach:" << endl << endl;
    file << "| \t\tWezel" << "\t\t\t\t\tWartosc dokladna" << "\t\tWartosc funkcji aproksymujacej" << endl;

    for (int i = 0; i < liczbawezlow; ++i)
    {
        file << "|\t\t";
        file << (*wezly)[i];
        file << "\t\t\t";
        file << (*fx)[i];
        file << "\t\t\t";
        file << (*qx)[i] << endl;
    }


    file.unsetf(ios::fixed);
    file.setf(ios::scientific);

    file << "| Wspolczynniki A" << endl;

    for (int i = 0; i <= stopien; ++i)
        file << "| " << i << ": " << (*A)[i] << endl;
    file << "+-------------------------------------------------------------------------------------------------------+" << endl << endl;
    file.close();
    return plik;
}

// Tworzy Excelowskie pliki dla kazdego przypadku stopnia
// Zwraca nazwę pliku
string _ToExcel(const string dir, const Wektor *wezly, const Wektor *fx, const Wektor *qx, const int liczbawezlow, const int stopien, const double granica_lewa, const double granica_prawa)
{
    ostringstream  ss_file;
    ss_file << "Stopien_" << stopien << "-liczbaWezlow_" << liczbawezlow << "-granicaL_" << granica_lewa << "-granicaP_" << granica_prawa << ".csv";
    string plik = ss_file.str();
    fstream file(dir + (string)DIR_DELIMITER + plik, ios::out | ios::trunc);

    // Precyzja
    file.precision(10);

    // Kropka na przecinki w liczbach double:
    file.imbue(locale(""));

    file << "\"x\"" << DELIMITER << "\"F(x)\"" << DELIMITER << "\"Q(x)\"" << endl;
    for (int i = 0; i < liczbawezlow; ++i)
        file << (*wezly)[i] << DELIMITER << (*fx)[i] << DELIMITER << (*qx)[i] << endl;

    file.close();
    return plik;

}

// Zapisuje posortowane bledy apro. wraz z całą resztą informacji
void _SaveLog(const string dir, const vector<DaneLog> logi)
{

    fstream file(dir + (string)DIR_DELIMITER + "bledy_apro_posortowane.txt", ios::out | ios::trunc);
    file.imbue(locale(""));
    file.precision(10);
    for (vector<DaneLog>::const_iterator it = logi.begin(); it != logi.end(); ++it)
        file << *it << endl;
    file.close();
}

int main()
{
    char breakout;

    // Pod Windowsem  koniecznie dwa: "D:" + "\\" + "katalog"
    string katalog_wynikowy = "D:" + (string)DIR_DELIMITER + "apro";

    // To można zmienić
    int liczbawezlow = 50;
    vector<int> stopienie = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 };
    vector<vector<double>> granice = { { -1, 1 }, { -2, 2 }, { -3, 3 }, { -4, 4 }, { -5, 5 }, { -10, 10 }, { -20, 20 } };

    //vector<int> stopienie = { 19 };
    //vector<vector<double>> granice = { { -1, 1 }, { -2, 2 }, { -3, 3 }, { -4, 4 }, { -5, 5 }, { -10, 10 }, { -20, 20 } };


    // Nie modyfikowac - tymczasowe zmienne:
    double granica_lewa;
    double granica_prawa;
    double blad;
    int stopien;
    string plik_raportu;
    string plik_excela;

    // Nasz Log
    vector<DaneLog> log;

    for (vector<vector<double>>::iterator it = granice.begin(); it != granice.end(); ++it)
    {
        vector<double> granice_tmp = *it;

         granica_lewa = granice_tmp[0];
         granica_prawa = granice_tmp[1];

        for each (stopien in stopienie)
        {
            blad = 0.0;
            Wektor *Fx = new Wektor();
            Wektor *Wezly = new Wektor();

            _Przygotuj(Fx, Wezly, liczbawezlow, granica_lewa, granica_prawa);

            Macierz *G = _G(Wezly, stopien, liczbawezlow);
            Wektor *R = _R(Wezly, Fx, stopien, liczbawezlow);

            Wektor *A = nullptr;
            Wektor *Qx = nullptr;


            if (!_MacierzTrj(R, G, stopien))
            {
                delete Fx, Wezly, G, R;
                cout << "Dzielenie przez zero przy stopniu: "<< stopien << ", granicy l.: "<< granica_lewa << " i p.: "<< granica_prawa << ". Wychodze..." << endl;
                cin >> breakout;
                exit(1);
            }

            A = _A(R, G, stopien);
            Qx = _Qx(A, Wezly, stopien, liczbawezlow);
            blad = _Blad(Qx, Fx, liczbawezlow);


            plik_raportu =_Raport(katalog_wynikowy, Wezly, Fx, Qx, A, liczbawezlow, stopien, granica_lewa, granica_prawa, blad);
            plik_excela = _ToExcel(katalog_wynikowy, Wezly, Fx, Qx, liczbawezlow, stopien, granica_lewa, granica_prawa);

            // Zapisujemy log
            DaneLog danelog(blad, granica_prawa, granica_lewa, stopien, plik_raportu, plik_excela);
            log.push_back(danelog);

            // Przygotowujemy się do następnej iteracji
            delete Fx, Wezly, G, R, A, Qx;

            Fx = NULL;
            Wezly = NULL;
            G = NULL;
            R = NULL;
            A = NULL;
            Qx = NULL;

        }
    }

    // sortujemy
    sort(log.begin(), log.end());

    // zapisujemy
    _SaveLog(katalog_wynikowy, log);


    return 0;
}