initial commit/wip

lukas-weber · Jan 9, 2023 · 200b877 · 200b877
commit 200b877
Show file tree

Hide file tree

Showing 9 changed files with 17,880 additions and 0 deletions.
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,3 @@
+build
+subprojects/*
+!subprojects/*.wrap
diff --git a/meson.build b/meson.build
@@ -0,0 +1,28 @@
+project('solp', 'cpp',
+        version : '0.1',
+        license : 'MIT')
+
+eigen_dep = dependency('eigen', fallback : ['eigen', 'eigen_dep'])
+
+subdir('src')
+
+solp_deps = [eigen_dep]
+
+if not meson.is_subproject()
+  pkg_mod = import('pkgconfig')
+  pkg_mod.generate(libraries : libsolp,
+                   version : meson.project_version(),
+                   name : 'lib'+meson.project_name(),
+                   filebase : meson.project_name(),
+                   description : 'A Simple Open-source Linear Programming library.')
+  install_headers(solp_headers)
+endif
+
+solp_dep = declare_dependency(
+  include_directories : [include_directories('src')],
+  link_with : libsolp,
+  dependencies : solp_deps,
+  version : meson.project_version()
+)
+
+subdir('test')
diff --git a/src/meson.build b/src/meson.build
@@ -0,0 +1,13 @@
+
+solp_srcs = files([
+  'solp.cpp',
+])
+
+solp_headers = files([
+  'solp.h',
+])
+
+libsolp = shared_library(meson.project_name(), solp_srcs,
+                         dependencies : [eigen_dep],
+                         install : not meson.is_subproject())
+
diff --git a/src/solp.cpp b/src/solp.cpp
@@ -0,0 +1,112 @@
+#include "solp.h"
+#include <iostream>
+#include <Eigen/Dense>
+
+namespace solp {
+
+result solve(const std::vector<double> &objective, const std::vector<constraint> &constraints) {
+	Eigen::VectorXd obj = Eigen::Map<const Eigen::VectorXd>(objective.data(), objective.size());
+
+	Eigen::MatrixXd A(constraints.size(), obj.size());
+	Eigen::VectorXd b;
+
+	for(int i = 0; i < A.rows(); i++) {
+		A.row(i) = Eigen::Map<const Eigen::RowVectorXd>(constraints[i].coeff.data(), constraints[i].coeff.size());
+		b(i) = constraints[i].rhs;
+	}
+
+	result r;
+	return r;
+
+}
+
+static std::vector<int> find_column_basis(Eigen::MatrixXd &A) {
+	auto &QR = A.colPivHouseholderQr();
+	auto &R = QR.matrixR().template triangularView<Eigen::Upper>();
+	auto &P = QR.colsPermutation();
+
+	assert(QR.rank() == A.rows());
+
+	std::vector<int> basis(R.rows());
+	for(int i = 0; i < R.rows(); i++) {
+		for(int j = i; j < R.cols(); j++) {
+			if(R(i,j) != 0) {
+				basis[i] = P.indices()(j);
+				break;
+			}
+		}
+	}
+
+	return basis;
+}
+
+static result canonical_simplex(Eigen::VectorXd &objective, Eigen::MatrixXd &A, Eigen::VectorXd &b) {
+	assert(A.cols() > A.rows());
+
+	int nb = A.rows();
+	int nn = A.cols()-A.rows();
+
+	Eigen::VectorXd xb(nb);
+	Eigen::VectorXd d(nb);
+
+	Eigen::VectorXd sn(nn);
+	Eigen::VectorXd lambda(nn);
+
+	//assert(A.leftCols(nb).isApprox(Eigen::MatrixXd::Identity(nb,nb))); // problem canonical?
+
+	std::vector<int> idxs(A.cols());
+	for(size_t i = 0; i < idxs.size(); i++) {
+		idxs[i] = i;
+	}
+
+
+	while(true) {
+		lambda = A.leftCols(nb).transpose().colPivHouseholderQr().solve(objective.head(nb));
+		sn = objective.tail(nn) - A.rightCols(nn).transpose()*lambda;
+
+		int q{-1};
+		int minidx = A.cols()+1;
+		for(int i = 0; i < nn; i++) {
+			if(sn(i) < 0 && idxs[nb+i] < minidx) {
+				q = nb + i;
+				minidx = idxs[q];
+			}
+		}
+
+		auto Bfac = A.leftCols(nb).colPivHouseholderQr();		
+		xb = Bfac.solve(b.head(nb)); 
+
+		if(q < 0) {
+			break;
+		}
+
+		d = Bfac.solve(A.col(q));
+
+		int p = -1;
+		double rmin = INFINITY; 
+		for(int i = 0; i < nb; i++) {
+			if(d(i) > 0 && xb(i)/d(i) < rmin) {
+				p = i;
+				rmin = xb(i)/d(i);
+			}
+		}
+		if(p < 0) {
+			throw std::runtime_error("problem unbounded");
+		}
+		std::cout << p << "<-" << q << "\n";
+
+		std::swap(idxs[p], idxs[q]);
+		A.col(p).swap(A.col(q));
+		std::swap(objective(p), objective(q));
+	}
+
+	result res;
+	res.x.resize(A.cols());
+	for(int i = 0; i < nb; i++) {
+		res.x[idxs[i]] = xb[i];
+		//res.basis[i] = idxs[i];	
+	}
+	return res;
+}
+
+}
diff --git a/src/solp.h b/src/solp.h
@@ -0,0 +1,40 @@
+#ifndef solp_h_INCLUDED
+#define solp_h_INCLUDED
+
+#include <vector>
+namespace solp {
+
+struct result {
+	std::vector<double> x;
+
+	// used internally
+	std::vector<int> basis;
+};
+
+// constraint represents an equality constraint on the variables x
+//
+//    coeff * x = rhs
+//    
+struct constraint {
+	std::vector<double> coeff;
+	double rhs;
+};
+
+// solve the linear programming problem
+//
+//     Minimize objective * x
+//     with
+//         x >= 0
+//         and the given constraints.
+//
+// Note that constraints only encode equality constraints. However,
+// every linear programming problem can be brought into this standard form
+// by shifting x and introducing additional slack variables.
+//
+// If there is no solution or the problem is unbounded, a solp::exception is thrown.
+// 
+result solve(const std::vector<double> &objective, const std::vector<constraint> &constraints);
+
+}
+#endif // solp_h_INCLUDED
+
diff --git a/subprojects/eigen.wrap b/subprojects/eigen.wrap
@@ -0,0 +1,10 @@
+[wrap-file]
+directory = eigen-eigen-323c052e1731
+source_url = http://bitbucket.org/eigen/eigen/get/3.3.7.tar.bz2
+source_filename = eigen-eigen-323c052e1731.tar.bz2
+source_hash = 9f13cf90dedbe3e52a19f43000d71fdf72e986beb9a5436dddcd61ff9d77a3ce
+
+patch_url = https://github.com/mesonbuild/eigen/releases/download/3.3.7-1/eigen.zip
+patch_filename = eigen-3.3.7-1-wrap.zip
+patch_hash = 3d8ffc134e8af95e5de6c5d4614971028d27bf5f01f3082179522d536750261b
+