#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <cmath>
#include <algorithm>
#include <stdlib.h>

using namespace std;
bool WczytajPlik(double ** &A, double * &b, unsigned *matrix_size);
void WyswietlMenu();
bool CzySymetryczna(double ** &A, double * &b, unsigned matrix_size);
void WyswietlMacierz(double ** A, double * b, unsigned matrix_size);
void KopiowanieTablicy(double ** &A, double * &b, double ** &C, double * &d, unsigned matrix_size);
void Choleski(double ** A, double * b, unsigned matrix_size);
void MacierzLY(double ** L, double * b, unsigned matrix_size, vector<double> &y);
void MacierzLTY(double ** LT, vector<double> &y, vector<double> &x, unsigned matrix_size);
void TranspozycjaMacierzyLY(double ** L, double ** LT, unsigned matrix_size);
void WyswietlMacierzBezWartosci(double ** A, unsigned matrix_size);

int main()
{
    double ** A;
    double * b;
    double ** C;
    double * d;
    int wybor;
    unsigned matrix_size;
    bool czyWczytany = false;
    for(;;)
    {
        WyswietlMenu();
        cout<<"Wybor => "; cin>>wybor;
        system("clear");
        switch(wybor)
        {
            case 1:
                if(WczytajPlik(C, d, &matrix_size)) {
                    czyWczytany = true;
                    cout<<"Plik zostal wczytany.\n";
                    WyswietlMacierz(C, d, matrix_size);
                    KopiowanieTablicy(A, b, C, d, matrix_size);
                }
                else
                    cout<<"Nie udalo sie wczytac pliku\n";
            break;
            case 2:
                if(czyWczytany) {
                    if(CzySymetryczna(A, b, matrix_size)) {
                        cout<<"Macierz jest symetryczna. Przechodze do obliczen..."<<endl;
                        Choleski(A, b, matrix_size);
                    }
                    else
                        cout<<"Macierz nie jest symetryczna.\nNalezy wczytac inny uklad."<<endl;
                    }
                else
                    cout<<"Nalezy wczytac plik przed rozpoczeciem programu";
            break;
            case 3:
                return 0;
            break;
            default:
                cout<<"Nalezy wybrac opcje z zakresu [1-3]\n";
            break;
        }
    }

    return 0;
}

void WyswietlMenu()
{
    cout<<"\n----MENU----\n";
    cout<<"1. Wczytanie danych z pliku.\n";
    cout<<"2. Metoda Choleskiego.\n";
    cout<<"3. Wyjscie z programu.\n";
    cout<<"\n";
}

bool WczytajPlik(double ** &A, double * &b, unsigned *matrix_size)
{
    ifstream source_file("test.csv");
    if (!source_file.is_open())
    {
        cout <<"The file has not been open!"<<endl;
        return false;
    }
    source_file >> *matrix_size;
    A = new double*[*matrix_size];
    A[0] = new double[*matrix_size**matrix_size];
    for(unsigned i = 1; i< *matrix_size; i++)
        A[i] = A[i-1] + *matrix_size;
    b = new double[*matrix_size];
    char semicolumn;
    for (unsigned i = 0; i < *matrix_size+1; i++)
        source_file >> semicolumn;

    for (unsigned i = 0; i < *matrix_size; i++)
    {
        for (unsigned j = 0; j < *matrix_size; j++)
        {
            source_file >> A[i][j];
            source_file >> semicolumn;
        }
        source_file >> semicolumn;
        source_file >> b[i];
    }
    source_file.close();
    return true;
}

void WyswietlMacierz(double ** A, double * b, unsigned matrix_size)
{
    cout<<endl;
    for(int i=0; i<matrix_size; i++) {
        for(int j=0; j<matrix_size; j++) {
            cout<<setw(10)<<A[i][j]<<" ";
        }
        cout<<"|"<<setw(10)<<b[i];
        cout<<"\n";
    }
    cout<<endl;

}
void WyswietlMacierzBezWartosci(double ** A, unsigned matrix_size)
{
    cout<<endl;
    for(int i=0; i<matrix_size; i++) {
        for(int j=0; j<matrix_size; j++) {
            cout<<setw(10)<<A[i][j]<<" ";
        }
        cout<<"\n";
    }
    cout<<endl;

}

void KopiowanieTablicy(double ** &A, double * &b, double ** &C, double * &d, unsigned matrix_size)
{
    b = new double[matrix_size];
    A = new double*[matrix_size];
    for(int i=0; i<matrix_size; i++)
        A[i] = new double[matrix_size];
    for(int i=0; i<matrix_size; i++)
        for(int j=0; j<matrix_size; j++)
            A[i][j] = C[i][j];
    for(int i=0; i<matrix_size; i++)
        b[i] = d[i];
}

bool CzySymetryczna(double ** &A, double * &b, unsigned matrix_size)
{
    for(int i=0; i<matrix_size; i++)
        for(int j=0; j<matrix_size; j++)
            if(A[i][j] != A[j][i])
                return false;
    return true;
}

void Choleski(double ** A, double * b, unsigned matrix_size)
{
    // s = kolumna
    double suma=0;
    vector<double> y; // Ly = b
    vector<double> x; // LT * x = y
    double **L = new double*[matrix_size]; //Macierz L
    for(int i=0; i<matrix_size; i++)
        L[i] = new double[matrix_size];
    double **LT = new double*[matrix_size]; //Macierz Transponowana L
    for(int i=0; i<matrix_size; i++)
        LT[i] = new double[matrix_size];

    cout<<"\n\n------------Macierz A-----------\n";
    WyswietlMacierz(A, b, matrix_size);

    cout<<"Macierz L  przed rozkladem.\n";
    WyswietlMacierz(L, b, matrix_size);
    for(int s=0; s<matrix_size; s++)
    {
        cout<<"----Kolumna nr"<<s+1<<"----"<<endl;
        for(int i=s; i<matrix_size; i++)
        {
            suma=0;
            cout<<"Wiersz = "<<i+1<<endl;
            if(s == i)
            {
                //cout<<"s == i iteracja: "<<i<<endl;
                for(int j=0; j<s; j++)
                {
                    suma+= pow(L[s][j], 2);
                    //cout<<"L["<<s<<"]["<<s<<"] = "<<L[s][j]<<endl;
                    //cout<<"Suma = "<<suma<<endl;
                }

                L[s][s] = sqrt(A[s][s] - suma);
                cout<<"Wstawiam "<<L[s][s]<<" w miejsce o wspolrzednych ("<<s+1<<","<<s+1<<")"<<endl;
                WyswietlMacierz(L, b, matrix_size);

            }
            else
            {
                //cout<<"else s==1 iteracja: "<<i<<endl;
                for(int j=0; j<s; j++) {
                    suma+= L[i][j]*L[s][j];
                    //cout<<"Suma = "<<suma<<endl;
                }
                L[i][s] = (A[i][s] - suma)/L[s][s] ;
                cout<<"Wstawiam "<<L[i][s]<<" w miejsce o wspolrzednych ("<<i+1<<","<<s+1<<")"<<endl;
                WyswietlMacierz(L, b, matrix_size);
            }
        }
    }
    cout<<"Macierz L po rozkladzie.\n";
    WyswietlMacierz(L, b, matrix_size);
    cout<<"\n--- Ly = b ---\n\n";
    MacierzLY(L, b, matrix_size, y);
    cout<<"\n--- LT * x = y ---\n";
    cout<<"\n1)Transponowanie macierzy L\n";
    TranspozycjaMacierzyLY(L, LT, matrix_size);
    cout<<"\n2)Wyliczam wartosci x\n";
    MacierzLTY(LT, y, x, matrix_size);
}

void MacierzLY(double ** L, double * b, unsigned matrix_size, vector<double> &y)
{
    double suma = 0;
    int pomocnyIterator = matrix_size-2;
    y.push_back(b[0] / L[0][0]);


    for(int i=1; i<matrix_size; i++)
    {
        suma = 0;
        for(int j=0; j<pomocnyIterator; j++)
        {
            suma += L[i][j]*y[j];
            //cout<<"Suma = "<<suma<<endl;
        }
        b[i] = b[i] - suma;
        y.push_back(b[i] / L[i][i]);
        pomocnyIterator++;
    }

    for(int i=0; i<matrix_size; i++)
        cout<<"y"<<i+1<<" = "<<y[i]<<endl;
}

void MacierzLTY(double ** LT, vector<double> &y, vector<double> &x, unsigned matrix_size)
{
    double suma = 0;
    int pomocnyIterator = matrix_size-2;
    int iteratorOdcietych = 0;
    int pom = matrix_size -1;
    WyswietlMacierzBezWartosci(LT, matrix_size);

    x.push_back( y[pom] / LT[pom][pom] );
    for(int i=matrix_size-2; i>=0; i--)
    {
        //cout<<"\n\nWiersz nr"<<i+1<<"\n\n";
        suma = 0;
        iteratorOdcietych = 0;
        for(int j=matrix_size-1; j>pomocnyIterator; j--)
        {
            //cout<<"Wchodze do petli for z iterarotem J\n";
            suma += LT[i][j]*x[iteratorOdcietych];
            iteratorOdcietych++;
            //cout<<"Suma = "<<suma<<endl;
        }
        y[i] = y[i] - suma;
        pomocnyIterator--;
        x.push_back( y[i] / LT[i][i] );
    }
    cout<<"\nRozwiazanie ukladu:\n\n";
    for(int i=0; i<matrix_size; i++)
        cout<<"x"<<i+1<<" = "<<x[i]<<endl;

}

void TranspozycjaMacierzyLY(double ** L, double ** LT, unsigned matrix_size)
{
    for(int i=0; i<matrix_size; i++)
        for(int j=0; j<matrix_size; j++)
            LT[i][j] = L[j][i];

    WyswietlMacierzBezWartosci(LT, matrix_size);
}