Orienteering/native/salesman.cpp

#include <cstdint>
#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{
	float** weights;
	uint8_t** adjacencyMatrix;
	float lowerBound;
	bool finished = false;
	std::vector<std::pair<int, int>> paths;
};

struct SplitState {
	bool wasSplittable = true;
	State firstState;
	State secondState;
};


double reduceWeights(float** 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(float** 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, uint8_t** 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) {
	uint8_t** adjacencyMatrixCopy = new uint8_t*[numberOfCities];
	float** weightsCopy = new float*[numberOfCities];

	for (int row = 0; row < numberOfCities; row++) {
		adjacencyMatrixCopy[row] = new uint8_t[numberOfCities];
		weightsCopy[row] = new float[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(uint8_t** 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 float*[numberOfCities];
	initialState.adjacencyMatrix = new uint8_t*[numberOfCities];
	for (int row = 0; row < numberOfCities; row++) {
		initialState.weights[row] = new float[numberOfCities];
		initialState.adjacencyMatrix[row] = new uint8_t[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;
	uint8_t** 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