Move general matrix to simplex and fix bugs

.vscode/settings.json

@@ -0,0 +1,48 @@
+{
+    "files.associations": {
+        "array": "cpp",
+        "atomic": "cpp",
+        "bit": "cpp",
+        "*.tcc": "cpp",
+        "cctype": "cpp",
+        "clocale": "cpp",
+        "cmath": "cpp",
+        "compare": "cpp",
+        "concepts": "cpp",
+        "cstdarg": "cpp",
+        "cstddef": "cpp",
+        "cstdint": "cpp",
+        "cstdio": "cpp",
+        "cstdlib": "cpp",
+        "cwchar": "cpp",
+        "cwctype": "cpp",
+        "deque": "cpp",
+        "string": "cpp",
+        "unordered_map": "cpp",
+        "vector": "cpp",
+        "exception": "cpp",
+        "algorithm": "cpp",
+        "functional": "cpp",
+        "iterator": "cpp",
+        "memory": "cpp",
+        "memory_resource": "cpp",
+        "numeric": "cpp",
+        "random": "cpp",
+        "string_view": "cpp",
+        "system_error": "cpp",
+        "tuple": "cpp",
+        "type_traits": "cpp",
+        "utility": "cpp",
+        "initializer_list": "cpp",
+        "iosfwd": "cpp",
+        "iostream": "cpp",
+        "istream": "cpp",
+        "limits": "cpp",
+        "new": "cpp",
+        "numbers": "cpp",
+        "ostream": "cpp",
+        "stdexcept": "cpp",
+        "streambuf": "cpp",
+        "typeinfo": "cpp"
+    }
+}

main.cpp

@@ -4,7 +4,8 @@
 #include "simplex.h"
 
 int main() {
-    Vector C = {5, 4, 0, 0};
+    // TODO: Initially should be positive
+    Vector C = {-5, -4, 0, 0, 0, 0};
     Matrix A = {{6, 4, 1, 0, 0, 0}, {1, 2, 0, 1, 0, 0}, {-1, 1, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 1}};
     Vector b = {24, 6, 1, 2};
 
@@ -12,11 +13,12 @@ int main() {
     if(result.state == bounded) {
 
         if(result.solution == nullptr) {
-            std::cout << "HUYN: Khong tim duoc nghiem" << std::endl;
+            // TODO: error
         }
         std::cout << *result.solution << std::endl;
         std::cout << result.objective_fucntion_value << std::endl;
 
+        delete result.solution;
     }
 
     return 0;

simplex.cpp

@@ -16,32 +16,35 @@ struct Result {
 };
 
 Result Simplex(Vector C, Matrix A, Vector b, double eps = 0.01, bool maximize=true) {
-
-    int m = A.getRows();
-
 	Result result{};
-	std::vector<int> basicVars(A.getColumns() - A.getRows());
+	Matrix generalMatrix = createGeneralMatrix(A, C, b);
+
+	std::vector<int> basicVars(generalMatrix.getRows());
 	basicVars[0] = -1;
 
-	for (int i = 1; i < basicVars.size(); i++) {
+	for (size_t i = 1; i < basicVars.size(); i++) {
 		basicVars[i] = static_cast<int>(basicVars.size()) + i;
 	}
 
 	while (true) {
 		//3
 		int pivot_column_index = 0;
-		pivot_column_index = max_index(C);
+		pivot_column_index = min_index(generalMatrix[0]);
 
-		if (A[0][pivot_column_index] >= 0) {
+		if (generalMatrix[0][pivot_column_index] >= 0) {
+			DestroyMatrix destroyedGeneralMatrix = destroyGeneralMatrix(generalMatrix);
+			Matrix _A = destroyedGeneralMatrix.A;
+			Vector _C = destroyedGeneralMatrix.C;
+			Vector _b = destroyedGeneralMatrix.b;
 
 			result.state = bounded;
 			result.solution = new Vector(C.size());
 			for (int i = 0; i < C.size(); i++) {
 				result.solution->operator[](i) = 0;
 			}
-			for (int i = 1; i < basicVars.size(); i++) {
+			for (size_t i = 1; i < basicVars.size(); i++) {
 				if (basicVars[i] <= C.size()) {
-					(*result.solution)[basicVars[i]] = b.size() - 1;
+					result.solution->operator[](basicVars[i]) = _b[i];
 				}
 			}
 			result.objective_function_value = b[0];
@@ -49,21 +52,22 @@ Result Simplex(Vector C, Matrix A, Vector b, double eps = 0.01, bool maximize=tr
 		}
 
 		//4
-		Vector ratio_vector(m);
+		Vector ratio_vector(generalMatrix.getRows());
-		for (int i = 0; i < m; i++) {
+		for (int i = 1; i < generalMatrix.getRows(); i++) {
-			ratio_vector[i] = b[i] / A[i][pivot_column_index];
+			if (generalMatrix[i][pivot_column_index] != 0) {
+            	ratio_vector[i] = generalMatrix[i][generalMatrix.getColumns() - 1] / generalMatrix[i][pivot_column_index];
+        	} else {
+            	ratio_vector[i] = 0;
+        	}
 		}
-		int pivot_row_index = min_index(ratio_vector);
+		ratio_vector[0] = 0;
+
+		int pivot_row_index = min_index_positive(ratio_vector);
 
 		basicVars[pivot_row_index] = pivot_column_index;
 
 		//5
-		FracturedMatrix fractured_matrix(elimination(A, C, b, pivot_column_index, pivot_row_index));
+		elimination(generalMatrix, pivot_row_index, pivot_column_index);
-		A = fractured_matrix.A;
-		C = fractured_matrix.C;
-		b = fractured_matrix.b;
-		pivot_column_index = fractured_matrix.pivot_column_index;
-		pivot_row_index = fractured_matrix.pivot_row_index;
 	}
 
 	return result;

tools/elimination.cpp

@@ -43,8 +43,9 @@ Matrix createGeneralMatrix(Matrix& A, Vector& C, Vector& b) {
         generalMatrix[0][i] = C[i];
     }
 
-    for (int j = 0; j < A.getRows(); j++) {
+    // For objective function (j=0) the value is set to zero automatically
-        generalMatrix[j][A.getColumns()] = b[j];
+    for (int j = 1; j < A.getRows() + 1; j++) {
+        generalMatrix[j][A.getColumns()] = b[j-1];
     }
 
     for (int i = 0; i < A.getRows(); i++) {
@@ -56,8 +57,7 @@ Matrix createGeneralMatrix(Matrix& A, Vector& C, Vector& b) {
     return generalMatrix;
 }
 
-FracturedMatrix elimination(Matrix A, Vector C, Vector b, int pivot_column_index, int pivot_row_index) {
+void elimination(Matrix& generalMatrix, int pivot_row_index, int pivot_column_index) {
-    Matrix generalMatrix = createGeneralMatrix(A, C, b);
 
     int rows = generalMatrix.getRows();
     int cols = generalMatrix.getColumns();
@@ -78,25 +78,5 @@ FracturedMatrix elimination(Matrix A, Vector C, Vector b, int pivot_column_index
             generalMatrix[i][j] -= pivotColumnCoefficient * generalMatrix[pivot_row_index][j];
         }
     }
-
-    DestroyMatrix destroyedMatrix = destroyGeneralMatrix(generalMatrix);
-
-    pivot_column_index = max_index(destroyedMatrix.C);
-    
-    Vector ratio_vector(A.getRows());
-
-    for (int i = 0; i < A.getRows(); i++) {
-        if (A[i][pivot_column_index] != 0) {
-            ratio_vector[i] = b[i] / A[i][pivot_column_index];
-        }
-        
-        else {
-            ratio_vector[i] = 0;
-        }
-    }
-
-    pivot_row_index = min_index(ratio_vector);
-
-    return {A, C, b, pivot_column_index, pivot_row_index};
 }
 

tools/elimination.h

@@ -23,7 +23,7 @@ struct DestroyMatrix {
 };
 
 DestroyMatrix destroyGeneralMatrix(Matrix& generalMatrix);
-FracturedMatrix elimination(Matrix A, Vector C, Vector b, int pivot_column_index, int pivot_row_index);
+void elimination(Matrix&, int pivot_row_index, int pivot_column_index);
 Matrix createGeneralMatrix(Matrix& A, Vector& C, Vector& b);
 
 

tools/math.cpp

@@ -20,6 +20,26 @@ double max(Vector vector) {
     return temp;
 }
 
+int min_index_positive(Vector vector) {
+    int i = 0;
+    while (i < vector.size() && vector[i] <= 0) {
+        ++i;
+    }
+    if (i >= vector.size()) {
+        throw std::runtime_error("No positive min found");
+    }
+    int temp_index = i;
+    double temp = vector[i];
+
+    for (int j = i + 1; j < vector.size(); j++) {
+        if (vector[j] < temp && vector[j] > 0) {
+            temp_index = j;
+            temp = vector[j];
+        }
+    }
+    return temp_index;
+}
+
 int min_index(Vector vector) {
     int temp_index = 0;
     double temp = vector[0];

tools/math.h

@@ -7,6 +7,7 @@
 double min(Vector vector);
 double max(Vector vector);
 int min_index(Vector vector);
+int min_index_positive(Vector vector);
 int max_index(Vector vector);
 double min(std::vector<double> array);
 double max(std::vector<double> array);