ProjetMAM/Pivot
Archived
4
1
Fork 0
Code Issues Pull Requests Actions Releases 2 Activity

solver rework + refactor

Browse Source
This commit is contained in:
Simon Pribylski
2024-03-06 21:24:15 +01:00
parent 259750a794
commit 99624d1b00
10 changed files with 191 additions and 163 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
src/Gauss.cpp
View File

@@ -34,12 +34,12 @@ static void SimplifyLine(Matrix& mat, std::size_t line, std::size_t pivot_line,
}
}
void GaussJordan(Matrix& mat, bool reduite, bool normalise) {
void GaussJordan(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise) {
int indice_ligne_pivot = -1;
for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {
for (std::size_t j = 0; j < a_Matrix.GetColumnCount(); j++) {
int indice_ligne_pivot_trouve = FirstNotNullElementIndexOnColumn(mat, j, indice_ligne_pivot + 1);
int indice_ligne_pivot_trouve = FirstNotNullElementIndexOnColumn(a_Matrix, j, indice_ligne_pivot + 1);
if (indice_ligne_pivot_trouve < 0) // colonne de 0
continue; // on regarde la prochaine colonne
@@ -47,20 +47,20 @@ void GaussJordan(Matrix& mat, bool reduite, bool normalise) {
indice_ligne_pivot++;
if (indice_ligne_pivot_trouve != indice_ligne_pivot) {
SwapLines(mat, indice_ligne_pivot_trouve, indice_ligne_pivot);
SwapLines(a_Matrix, indice_ligne_pivot_trouve, indice_ligne_pivot);
}
Matrix::Element pivot = mat.at(indice_ligne_pivot, j);
Matrix::Element pivot = a_Matrix.at(indice_ligne_pivot, j);
if (normalise) {
DivideLine(mat, indice_ligne_pivot, pivot);
if (a_Normalise) {
DivideLine(a_Matrix, indice_ligne_pivot, pivot);
}
// On simplifie les autres lignes
for (std::size_t i = (reduite ? 0 : j); i < mat.GetRawCount(); i++) {
for (std::size_t i = (a_Reduite ? 0 : j); i < a_Matrix.GetRawCount(); i++) {
// Pour les lignes autre que la ligne pivot
if (i != static_cast<std::size_t>(indice_ligne_pivot)) {
SimplifyLine(mat, i, indice_ligne_pivot, j);
SimplifyLine(a_Matrix, i, indice_ligne_pivot, j);
}
}
}

87
src/Matrix.cpp
View File

@@ -7,29 +7,29 @@
#include <fstream>
#include <iostream>
Matrix::Matrix(std::size_t lignes, std::size_t colonnes) : m_Raws(lignes), m_Columns(colonnes) {
Matrix::Matrix(std::size_t a_Raws, std::size_t a_Columns) : m_Raws(a_Raws), m_Columns(a_Columns) {
m_Data.resize(m_Raws * m_Columns);
}
Matrix::Matrix(std::size_t lignes, std::size_t colonnes, std::initializer_list<Element>&& initList) :
m_Raws(lignes), m_Columns(colonnes) {
m_Data = initList;
Matrix::Matrix(std::size_t a_Raws, std::size_t a_Columns, std::initializer_list<Element>&& a_Elements) :
m_Raws(a_Raws), m_Columns(a_Columns) {
m_Data = a_Elements;
m_Data.resize(m_Raws * m_Columns);
}
Matrix Matrix::operator*(const Matrix& other) const {
if (m_Columns != other.m_Raws) {
Matrix Matrix::operator*(const Matrix& a_Other) const {
if (m_Columns != a_Other.m_Raws) {
std::cerr << "Mutiplication impossible car la dimensions des matrices est incompatible" << std::endl;
return {};
}
Matrix result(m_Raws, other.m_Columns);
Matrix result(m_Raws, a_Other.m_Columns);
for (std::size_t i = 0; i < m_Raws; ++i) {
for (std::size_t j = 0; j < other.m_Columns; ++j) {
for (std::size_t j = 0; j < a_Other.m_Columns; ++j) {
Element sum = 0;
for (std::size_t k = 0; k < m_Columns; k++) {
sum += at(i, k) * other.at(k, j);
sum += at(i, k) * a_Other.at(k, j);
}
result.at(i, j) = sum;
}
@@ -47,35 +47,35 @@ void Matrix::Transpose() {
*this = result;
}
Matrix Matrix::Identity(std::size_t taille) {
Matrix id {taille, taille};
for (std::size_t i = 0; i < taille; i++) {
for (std::size_t j = i; j < taille; j++) {
Matrix Matrix::Identity(std::size_t a_Size) {
Matrix id {a_Size, a_Size};
for (std::size_t i = 0; i < a_Size; i++) {
for (std::size_t j = i; j < a_Size; j++) {
id.at(i, j) = (i == j);
}
}
return id;
}
Matrix Matrix::ColumnVector(std::initializer_list<Element>&& initList) {
Matrix result {initList.size(), 1};
Matrix Matrix::ColumnVector(std::initializer_list<Element>&& a_Elements) {
Matrix result {a_Elements.size(), 1};
result.m_Data = initList;
result.m_Data = a_Elements;
return result;
}
Matrix Matrix::RawVector(std::initializer_list<Element>&& initList) {
Matrix result {1, initList.size()};
Matrix Matrix::RawVector(std::initializer_list<Element>&& a_Elements) {
Matrix result {1, a_Elements.size()};
result.m_Data = initList;
result.m_Data = a_Elements;
return result;
}
void Matrix::Augment(const Matrix& droite) {
assert(droite.m_Raws == m_Raws);
Matrix temp {m_Raws, m_Columns + droite.m_Columns};
void Matrix::Augment(const Matrix& a_Right) {
assert(a_Right.m_Raws == m_Raws);
Matrix temp {m_Raws, m_Columns + a_Right.m_Columns};
for (std::size_t i = 0; i < m_Raws; i++) {
for (std::size_t j = 0; j < m_Columns; j++) {
@@ -84,49 +84,49 @@ void Matrix::Augment(const Matrix& droite) {
}
for (std::size_t i = 0; i < m_Raws; i++) {
for (std::size_t j = 0; j < droite.m_Columns; j++) {
temp.at(i, j + m_Columns) = droite.at(i, j);
for (std::size_t j = 0; j < a_Right.m_Columns; j++) {
temp.at(i, j + m_Columns) = a_Right.at(i, j);
}
}
*this = temp;
}
Matrix Matrix::operator+(const Matrix& other) const {
assert(GetColumnCount() == other.GetColumnCount() && GetRawCount() == other.GetRawCount());
Matrix Matrix::operator+(const Matrix& a_Other) const {
assert(GetColumnCount() == a_Other.GetColumnCount() && GetRawCount() == a_Other.GetRawCount());
Matrix result = *this;
for (std::size_t i = 0; i < GetRawCount(); i++) {
for (std::size_t j = 0; j < GetColumnCount(); j++) {
result.at(i, j) += other.at(i, j);
result.at(i, j) += a_Other.at(i, j);
}
}
return result;
}
Matrix Matrix::operator-(const Matrix& other) const {
assert(GetColumnCount() == other.GetColumnCount() && GetRawCount() == other.GetRawCount());
Matrix Matrix::operator-(const Matrix& a_Other) const {
assert(GetColumnCount() == a_Other.GetColumnCount() && GetRawCount() == a_Other.GetRawCount());
Matrix result = *this;
for (std::size_t i = 0; i < GetRawCount(); i++) {
for (std::size_t j = 0; j < GetColumnCount(); j++) {
result.at(i, j) -= other.at(i, j);
result.at(i, j) -= a_Other.at(i, j);
}
}
return result;
}
bool Matrix::operator==(const Matrix& other) const {
if (m_Raws != other.m_Raws || m_Columns != other.m_Columns)
bool Matrix::operator==(const Matrix& a_Other) const {
if (m_Raws != a_Other.m_Raws || m_Columns != a_Other.m_Columns)
return false;
for (std::size_t i = 0; i < m_Raws; i++) {
for (std::size_t j = 0; j < m_Columns; j++) {
if (!IsEqualZero(at(i, j) - other.at(i, j)))
if (!IsEqualZero(at(i, j) - a_Other.at(i, j)))
return false;
}
}
@@ -134,12 +134,12 @@ bool Matrix::operator==(const Matrix& other) const {
return true;
}
Matrix::Element& Matrix::at(std::size_t ligne, std::size_t colonne) {
return m_Data[ligne * m_Columns + colonne];
Matrix::Element& Matrix::at(std::size_t a_Raw, std::size_t a_Column) {
return m_Data[a_Raw * m_Columns + a_Column];
}
Matrix::Element Matrix::at(std::size_t ligne, std::size_t colonne) const {
return m_Data[ligne * m_Columns + colonne];
Matrix::Element Matrix::at(std::size_t a_Raw, std::size_t a_Column) const {
return m_Data[a_Raw * m_Columns + a_Column];
}
std::size_t Matrix::GetRawCount() const {
@@ -150,14 +150,15 @@ std::size_t Matrix::GetColumnCount() const {
return m_Columns;
}
Matrix Matrix::SubMatrix(std::size_t origine_ligne, std::size_t origine_colonne, std::size_t ligne, std::size_t colonne) const {
assert(m_Raws >= origine_ligne + ligne && m_Columns >= origine_colonne + colonne);
Matrix Matrix::SubMatrix(
std::size_t a_RawOrigin, std::size_t a_ColumnOrigin, std::size_t a_RawCount, std::size_t a_ColumnCount) const {
assert(m_Raws >= a_RawOrigin + a_RawCount && m_Columns >= a_ColumnOrigin + a_ColumnCount);
Matrix result {ligne, colonne};
Matrix result {a_RawCount, a_ColumnCount};
for (std::size_t i = 0; i < ligne; i++) {
for (std::size_t j = 0; j < colonne; j++) {
result.at(i, j) = at(i + origine_ligne, j + origine_colonne);
for (std::size_t i = 0; i < result.GetRawCount(); i++) {
for (std::size_t j = 0; j < result.GetColumnCount(); j++) {
result.at(i, j) = at(i + a_RawOrigin, j + a_ColumnOrigin);
}
}

27
src/Solver.cpp
View File

@@ -2,10 +2,8 @@
#include "Gauss.h"
Solver::Solver(const Matrix& mat) : m_Matrix(mat) {}
Vect Solver::Image() const {
Matrix result = m_Matrix;
Vect Solver::Image(const Matrix& a_Matrix) const {
Matrix result = a_Matrix;
result.Transpose();
Gauss::GaussJordan(result, true, true);
result.Transpose();
@@ -13,27 +11,28 @@ Vect Solver::Image() const {
}
// https://en.wikipedia.org/wiki/Kernel_(linear_algebra)#Computation_by_Gaussian_elimination
Vect Solver::Kernel() const {
Matrix result = m_Matrix;
Vect Solver::Kernel(const Matrix& a_Matrix) const {
Matrix result = a_Matrix;
result.Transpose();
result.Augment(Matrix::Identity(result.GetRawCount()));
Gauss::GaussJordan(result, true, true);
result.Transpose();
// nombre de colonnes non nulles
std::size_t origine_colonne = Vect(result.SubMatrix(0, 0, m_Matrix.GetRawCount(), m_Matrix.GetColumnCount())).GetCardinal();
std::size_t origine_colonne = Vect(result.SubMatrix(0, 0, a_Matrix.GetRawCount(), a_Matrix.GetColumnCount())).GetCardinal();
return {result.SubMatrix(m_Matrix.GetRawCount(), origine_colonne, result.GetRawCount() - m_Matrix.GetRawCount(),
return {result.SubMatrix(a_Matrix.GetRawCount(), origine_colonne, result.GetRawCount() - a_Matrix.GetRawCount(),
result.GetColumnCount() - origine_colonne)};
}
VectAffine Solver::TriangularSystem() const {
Matrix mat = m_Matrix;
VectAffine Solver::RectangularSystem(const Matrix& a_MatrixA, const Matrix& a_VectorB) const {
Matrix mat = a_MatrixA;
mat.Augment(a_VectorB);
Gauss::GaussJordan(mat, true, true);
Solver solver {mat.SubMatrix(0, 0, mat.GetRawCount(), mat.GetColumnCount() - 1)};
Solver solver;
Vect noyau = solver.Kernel();
Vect noyau = solver.Kernel(a_MatrixA);
Matrix origin = mat.SubMatrix(0, mat.GetColumnCount() - 1, mat.GetRawCount(), 1);
// on rajoute des 0 si il faut
@@ -50,6 +49,6 @@ VectAffine Solver::TriangularSystem() const {
return {noyau, fullOrigin};
}
std::size_t Solver::Rank() const {
return Image().GetCardinal();
std::size_t Solver::Rank(const Matrix& a_Matrix) const {
return Image(a_Matrix).GetCardinal();
}

44
src/Vect.cpp
View File

@@ -14,7 +14,7 @@ static bool IsColumnNull(Matrix& mat, std::size_t column) {
return true;
}
Vect::Vect(const Matrix& mat) : m_Data(mat) {
Vect::Vect(const Matrix& a_Matrix) : m_Data(a_Matrix) {
Simplify();
}
@@ -29,50 +29,50 @@ void Vect::Simplify() {
m_Data = mat;
}
Matrix Vect::GetVector(std::size_t index) const {
return m_Data.SubMatrix(0, index, m_Data.GetRawCount(), 1);
Matrix Vect::GetVector(std::size_t a_Index) const {
return m_Data.SubMatrix(0, a_Index, m_Data.GetRawCount(), 1);
}
std::size_t Vect::GetCardinal() const {
return m_Data.GetColumnCount();
}
bool Vect::IsElementOf(const Matrix& vec) const {
bool Vect::IsElementOf(const Matrix& a_Vector) const {
Vect base = *this;
base.AddVector(vec);
base.AddVector(a_Vector);
return base.GetCardinal() == GetCardinal();
}
bool Vect::operator==(const Vect& other) const {
if (GetDimension() != other.GetDimension() || GetCardinal() != other.GetCardinal())
bool Vect::operator==(const Vect& a_Other) const {
if (GetDimension() != a_Other.GetDimension() || GetCardinal() != a_Other.GetCardinal())
return false;
// on vérifie si chaque vecteur de la deuxième base appartient à l'espace vectoriel engendré par la première base
for (std::size_t i = 0; i < GetCardinal(); i++) {
if (!IsElementOf(other.GetVector(i)))
if (!IsElementOf(a_Other.GetVector(i)))
return false;
}
return true;
}
void Vect::AddVector(const Matrix& mat) {
m_Data.Augment(mat);
void Vect::AddVector(const Matrix& a_Vector) {
m_Data.Augment(a_Vector);
m_Data.Transpose();
Gauss::GaussJordan(m_Data, false, false);
m_Data.Transpose();
Simplify();
}
bool Vect::operator!=(const Vect& other) const {
return !(*this == other);
bool Vect::operator!=(const Vect& a_Other) const {
return !(*this == a_Other);
}
Matrix Vect::GetLinearSystem() const {
Matrix vect = m_Data;
vect.Transpose();
Solver solver {vect};
vect = solver.Kernel().m_Data;
Solver solver;
vect = solver.Kernel(vect).m_Data;
vect.Transpose();
return vect;
}
@@ -81,9 +81,17 @@ std::size_t Vect::GetDimension() const {
return m_Data.GetRawCount();
}
VectAffine::VectAffine(const Vect& base, const Matrix& origine) :
m_Base(base), m_Origin(origine.SubMatrix(0, 0, m_Base.GetDimension(), 1)) {}
VectAffine::VectAffine(const Vect& a_Base, const Matrix& a_Origin) :
m_Base(a_Base), m_Origin(a_Origin.SubMatrix(0, 0, m_Base.GetDimension(), 1)) {}
bool VectAffine::IsElementOf(const Matrix& vec) const {
return m_Base.IsElementOf(vec - m_Origin);
bool VectAffine::IsElementOf(const Matrix& a_Vector) const {
return m_Base.IsElementOf(a_Vector - m_Origin);
}
Matrix VectAffine::GetLinearSystem() const {
Matrix result = m_Base.GetLinearSystem();
result.Augment(m_Origin.SubMatrix(0, 0, result.GetRawCount(), 1));
return result;
}

8
src/main.cpp
View File

@@ -23,10 +23,10 @@ void test() {
Matrix mat2 = LoadMatrix("matrice4x4.mat");
Print(mat2);
Solver solver {mat2};
Solver solver;
Vect image = solver.Image();
Vect noyau = solver.Kernel();
Vect image = solver.Image(mat2);
Vect noyau = solver.Kernel(mat2);
std::cout << "\tImage :\n";
Print(image);
@@ -38,7 +38,7 @@ void test() {
Print(noyau.GetLinearSystem());
std::cout << "\n\n";
Print(solver.TriangularSystem());
// Print(solver.TriangularSystem(mat2));
}
void prompt() {