initial commit

2026-08-29 09:29:39 +02:00
commit ccbf9c12a0
16 changed files with 3416 additions and 0 deletions
@@ -0,0 +1,425 @@
+#ifdef __EMSCRIPTEN__
+#include <emscripten/bind.h>
+#include <emscripten/val.h>
+#endif
+
+#include <cstdio>
+#include <ctime>
+#include <queue>
+#include <utility>
+#include <vector>
+
+struct City {
+	double x;
+	double y;
+};
+
+struct State{
+	double** weights;
+	int** adjacencyMatrix;
+	double lowerBound;
+	bool finished = false;
+	std::vector<std::pair<int, int>> paths;
+};
+
+struct SplitState {
+	bool wasSplittable = true;
+	State firstState;
+	State secondState;
+};
+
+
+double reduceWeights(double** weights, int numberOfCities) {
+	double totalReduced = 0;
+
+	for (int row = 0; row < numberOfCities; row++) {
+		double smallestOnRow = 1e100;
+		for (int column = 0; column < numberOfCities; column++) {
+			if (weights[row][column] < smallestOnRow) {
+				smallestOnRow = weights[row][column];
+			}
+		}
+
+		if (smallestOnRow > 1e50) {
+			continue;
+		}
+
+		totalReduced += smallestOnRow;
+
+		for (int column = 0; column < numberOfCities; column++) {
+			weights[row][column] -= smallestOnRow;
+		}
+	}
+
+	for (int column = 0; column < numberOfCities; column++) {
+		double smallestOnColumn = 1e100;
+		for (int row = 0; row < numberOfCities; row++) {
+			if (weights[row][column] < smallestOnColumn) {
+				smallestOnColumn = weights[row][column];
+			}
+		}
+
+		if (smallestOnColumn > 1e50) {
+			continue;
+		}
+		
+		totalReduced += smallestOnColumn;
+
+		for (int row = 0; row < numberOfCities; row++) {
+			weights[row][column] -= smallestOnColumn;
+		}
+	}
+
+	return totalReduced;
+}
+
+std::pair<int, int> findPivotPoint(double** weights, int numberOfCities) {
+	int bestRow = -1;
+	int bestColumn = -1;
+	double bestIncrease = -1;
+
+	for (int testRow = 0; testRow < numberOfCities; testRow++) {
+		for (int testColumn = 0; testColumn < numberOfCities; testColumn++) {
+			// Only look for pivot cells that have a value of zero
+			if (weights[testRow][testColumn] > 0.0001) {
+				continue;
+			}
+
+			double smallestOnRow = 1e100;
+			double smallestOnColumn = 1e100;
+
+			for (int otherColumn = 0; otherColumn < numberOfCities; otherColumn++) {
+				if (otherColumn == testColumn) {
+					continue;
+				}
+
+				if (weights[testRow][otherColumn] < smallestOnRow) {
+					smallestOnRow = weights[testRow][otherColumn];
+				}
+			}
+
+			for (int otherRow = 0; otherRow < numberOfCities; otherRow++) {
+				if (otherRow == testRow) {
+					continue;
+				}
+
+				if (weights[otherRow][testColumn] < smallestOnColumn) {
+					smallestOnColumn = weights[otherRow][testColumn];
+				}
+			}
+
+			double totalIncrease = smallestOnRow + smallestOnColumn;
+			if (totalIncrease > bestIncrease) {
+				bestIncrease = totalIncrease;
+				bestRow = testRow;
+				bestColumn = testColumn;
+			}
+		}
+	}
+
+	return std::pair<int, int>(bestRow, bestColumn);
+}
+
+int getDegreeOfNode(int nodeNumber, int** adjacencyMatrix, int numberOfCities) {
+	int degree = 0;
+	for (int column = 0; column < numberOfCities; column++) {
+		degree += adjacencyMatrix[nodeNumber][column];
+	}
+	return degree;
+}
+
+void multiplyAdjacencyMatrices(int** matrix, int** multiplier, int size) {
+	int** temporary = new int*[size];
+	for (int row = 0; row < size; row++) {
+		temporary[row] = new int[size];
+		for (int column = 0; column < size; column++) {
+			int result = 0;
+			for (int k = 0; k < size; k++) {
+				result += matrix[row][k] * multiplier[k][column];
+			}
+			temporary[row][column] = result;
+		}
+	}
+
+	// Copy into original matrix
+	for (int row = 0; row < size; row++) {
+		for (int column = 0; column < size; column++) {
+			matrix[row][column] = temporary[row][column];
+		}
+		delete [] temporary[row];
+	}
+	delete [] temporary;
+}
+
+void disallowSubloops(State &state, int pivotRow, int pivotColumn, int numberOfCities) {
+	// All points with degree one are endpoint
+	int totalDegree = 0;
+	for (int node = 0; node < numberOfCities; node++) {
+		int degree = getDegreeOfNode(node, state.adjacencyMatrix, numberOfCities);
+		totalDegree += degree;
+		if (degree == 2) {
+			for (int otherNode = 0; otherNode < numberOfCities; otherNode++) {
+				state.weights[node][otherNode] = 1e100;
+				state.weights[otherNode][node] = 1e100;
+			}
+		}
+	}
+
+	// Add this path to the existing paths
+	std::pair<int, int> path(pivotRow, pivotColumn);
+	auto iterator = state.paths.begin();
+	while (iterator != state.paths.end()) {
+		std::pair<int, int> otherPath = *iterator;
+
+		bool expandsPath = false;
+		if (path.first == otherPath.first) {
+			expandsPath = true;
+			path = std::pair<int,int>(path.second, otherPath.second);
+		} else if (path.first == otherPath.second) {
+			expandsPath = true;
+			path = std::pair<int,int>(path.second, otherPath.first);
+		} else if (path.second == otherPath.first) {
+			expandsPath = true;
+			path = std::pair<int,int>(path.first, otherPath.second);
+		} else if (path.second == otherPath.second) {
+			expandsPath = true;
+			path = std::pair<int,int>(path.first, otherPath.first);
+		}
+
+		if (expandsPath) {
+			iterator = state.paths.erase(iterator);
+		} else {
+			iterator++;
+		}
+	}
+
+	state.paths.push_back(path);
+
+	// There needs to be n edges in a loop, and so the total degree should be 2*n. If we're getting close to this, don't
+	// block of the ability to finish a loop
+	if (state.paths.size() == 1 && totalDegree == 2*(numberOfCities - 1)) {
+		// Finish the loop
+		int startingNode = state.paths.at(0).first;
+		int endingNode = state.paths.at(0).second;
+
+		state.adjacencyMatrix[startingNode][endingNode] = 1;
+		state.adjacencyMatrix[endingNode][startingNode] = 1;
+
+		state.finished = true;
+
+		return;
+	}
+
+	for (std::pair<int, int> path : state.paths) {
+		state.weights[path.first][path.second] = 1e100;
+		state.weights[path.second][path.first] = 1e100;
+	}
+}
+
+State createNewStateOnPivot(int pivotRow, int pivotColumn, State originalState, int numberOfCities) {
+	int** adjacencyMatrixCopy = new int*[numberOfCities];
+	double** weightsCopy = new double*[numberOfCities];
+
+	for (int row = 0; row < numberOfCities; row++) {
+		adjacencyMatrixCopy[row] = new int[numberOfCities];
+		weightsCopy[row] = new double[numberOfCities];
+		for (int column = 0; column < numberOfCities; column++) {
+			adjacencyMatrixCopy[row][column] = originalState.adjacencyMatrix[row][column];
+			weightsCopy[row][column] = originalState.weights[row][column];
+		}
+	}
+
+	adjacencyMatrixCopy[pivotRow][pivotColumn] = 1;
+	adjacencyMatrixCopy[pivotColumn][pivotRow] = 1;
+
+	for (int column = 0; column < numberOfCities; column++) {
+		weightsCopy[pivotRow][column] = 1e100;
+	}
+
+	for (int row = 0; row < numberOfCities; row++) {
+		weightsCopy[row][pivotColumn] = 1e100;
+	}
+
+	weightsCopy[pivotColumn][pivotRow] = 1e100;
+
+	State pivotedState;
+	pivotedState.adjacencyMatrix = adjacencyMatrixCopy;
+	pivotedState.weights = weightsCopy;
+	pivotedState.lowerBound = originalState.lowerBound;
+	for (std::pair<int, int> path : originalState.paths) {
+		pivotedState.paths.push_back(std::pair<int,int>(path.first, path.second));
+	}
+
+	disallowSubloops(pivotedState, pivotRow, pivotColumn, numberOfCities);
+
+	return pivotedState;
+}
+
+SplitState splitState(State state, int numberOfCities) {
+	SplitState splitState;
+
+	double reduction = reduceWeights(state.weights, numberOfCities);
+	state.lowerBound += reduction;
+
+	std::pair<int, int> pivotPoint = findPivotPoint(state.weights, numberOfCities);
+	if (pivotPoint.first == -1 || pivotPoint.second == -1) {
+		splitState.wasSplittable = false;
+
+		for (int row  = 0; row < numberOfCities; row++) {
+			delete [] state.adjacencyMatrix[row];
+			delete [] state.weights[row];
+		}
+		
+		delete [] state.adjacencyMatrix;
+		delete [] state.weights;
+
+		return splitState;
+	}
+	State pivotedState = createNewStateOnPivot(pivotPoint.first, pivotPoint.second, state, numberOfCities);
+
+	// Disallow the pivot in the original state
+	state.weights[pivotPoint.first][pivotPoint.second] = 1e100;
+
+	splitState.firstState = pivotedState;
+	splitState.secondState = state;
+	return splitState;
+}
+
+double getActualWeight(int** adjacencyMatrix, double** weights, int numberOfCities) {
+	double totalWeight = 0;
+	for (int row = 0; row < numberOfCities; row++) {
+		for (int column = 0; column < numberOfCities; column++) {
+			totalWeight += adjacencyMatrix[row][column] * weights[row][column];
+		}
+	}
+
+	return totalWeight / 2;
+}
+
+class StateComparator {
+	public:
+		bool operator()(State one, State two) {
+			return one.lowerBound > two.lowerBound;
+		}
+};
+
+std::vector<int> findShortestPath(std::vector<std::vector<double>> jsWeights) {
+
+	int numberOfCities = jsWeights.size();
+
+	double** weights = new double*[numberOfCities];
+	for (int row = 0; row < numberOfCities; row++) {
+		weights[row] = new double[numberOfCities];
+		for (int column = 0; column < numberOfCities; column++) {
+			weights[row][column] = jsWeights.at(row).at(column);
+			if (column == row) {
+				weights[row][column] = 1e100;
+			}
+		}
+	}
+	
+
+	// Create a state
+	State initialState;
+	initialState.weights = new double*[numberOfCities];
+	initialState.adjacencyMatrix = new int*[numberOfCities];
+	for (int row = 0; row < numberOfCities; row++) {
+		initialState.weights[row] = new double[numberOfCities];
+		initialState.adjacencyMatrix[row] = new int[numberOfCities];
+		for (int column = 0; column < numberOfCities; column++) {
+			initialState.adjacencyMatrix[row][column] = 0;
+			initialState.weights[row][column] = weights[row][column];
+		}
+	}
+	initialState.lowerBound = 0;
+
+	std::priority_queue<State, std::vector<State>, StateComparator> queue;
+	queue.push(initialState);
+
+	double currentBest = 1e100;
+	int** bestAdjacenyMatrix = NULL;
+
+	while (!queue.empty()) {
+		State nextState = queue.top();
+		queue.pop();
+
+		if (nextState.lowerBound > currentBest) {
+			break;
+		}
+
+		if (nextState.finished) {
+			double actualWeight = getActualWeight(nextState.adjacencyMatrix, weights, numberOfCities);
+			if (actualWeight < currentBest) {
+				currentBest = actualWeight;
+				bestAdjacenyMatrix = nextState.adjacencyMatrix;
+			} else {
+				// Delete and clean up memory
+				for (int row = 0; row < numberOfCities; row++) {
+					delete [] nextState.adjacencyMatrix[row];
+				}
+				delete [] nextState.adjacencyMatrix;
+			}
+
+			for (int row = 0; row < numberOfCities; row++) {
+				delete [] nextState.weights[row];
+			}
+			delete [] nextState.weights;
+
+			continue;
+		}
+
+		// If we're not finished, split the state and add the new ones to the queue
+		SplitState split = splitState(nextState, numberOfCities);
+		if (split.wasSplittable) {
+			queue.push(split.firstState);
+			queue.push(split.secondState);
+		}
+	}
+
+	// Clear up remaining queue for memory
+	while (!queue.empty()) {
+		State stateToDelete = queue.top();
+		queue.pop();
+
+		for (int row = 0; row < numberOfCities; row++) {
+			delete [] stateToDelete.adjacencyMatrix[row];
+			delete [] stateToDelete.weights[row];
+		}
+		delete [] stateToDelete.adjacencyMatrix;
+		delete [] stateToDelete.weights;
+	}
+
+	std::vector<int> path;
+
+	int current = 0;
+	int previous = -1;
+
+	path.push_back(current);
+
+	bool firstTime = true;
+	while (current != 0 || firstTime) {
+		firstTime = false;
+
+		for (int column = 0; column < numberOfCities; column++) {
+			if (bestAdjacenyMatrix[current][column] == 1 && column != previous) {
+				previous = current;
+				current = column;
+				break;
+			}
+		}
+
+		path.push_back(current);
+	}
+
+	return path;
+}
+
+#ifdef __EMSCRIPTEN__
+EMSCRIPTEN_BINDINGS(my_module) {
+	emscripten::register_vector<double>("WeightsRow");
+	emscripten::register_vector<std::vector<double>>("Weights");
+	emscripten::register_vector<int>("Path");
+
+	emscripten::function("findShortestPath", &findShortestPath);
+}
+#endif