GPU621/DPS921 T-eaSyPeasy

GPU621/DPS921 | Participants | Groups and Projects | Resources | Glossary

T-eaSyPeasy

Team Members

1. Berwout de Vries Robles

Introduction

For my project instead of picking a new method to learn to parallelise with, I decided to try to paralellise an algorithm that is fairly hard to write, to see if after following this course I could actually parallelise an algorithm. The thing is, in the workshops we have only been parallelising algorithms that have been specifically written to be parallelised. I wanted to change that and see if I could still do it. The algorithm I am parallelising is a Branch & Bound Algorithm by Little from 1966. A good explanation of the algorithm can be found here: https://www.youtube.com/watch?v=nN4K8xA8ShM, the explanation of the branch and bound I wrote starts at 34:30.

Source

#include <stdint.h>
#include <vector>
#include <algorithm>
#include <fstream>
#include <iostream>
#include <cmath>
#include <time.h>


struct Selection {
    uint16_t row;
    uint16_t col;
    uint16_t val;
    Selection(uint16_t x, uint16_t y, uint16_t v);

    bool operator==(const Selection& rhs)
    {
        return (col == rhs.col && row == rhs.row);
    }
    bool operator!=(const Selection& rhs)
    {
        return !(col == rhs.col && row == rhs.row);
    }
    Selection() : row(UINT16_MAX), col(UINT16_MAX), val(UINT16_MAX) {}
};

Selection::Selection(uint16_t r, uint16_t c, uint16_t v) {
    row = r;
    col = c;
    val = v;
}

struct Zero {
    uint16_t row;
    uint16_t col;
    uint16_t horPen;
    uint16_t vertPen;
    Zero(uint16_t row, uint16_t col);
    Zero() {};
};

Zero::Zero(uint16_t r, uint16_t c) {
    row = r;
    col = c;
};

struct Node {
    int lowerBound;
    int matrixSize;
    bool solutionFound = false;
    bool invalidSolution = false;
    uint16_t* originalMatrix;
    uint16_t* workingCopy;
    Node(int amountOfCities, uint16_t* matrix, std::vector<std::vector<Selection>> solution, uint16_t lb);
    Node(int amountOfCities, uint16_t* matrix, std::vector<std::vector<Selection>> solution);
    std::vector<std::vector<Selection>> solution;
    Node* left;
    Node* right;

    void evaluate(){
        std::vector<Selection> rowMinima;
        std::vector<Selection> colMinima;
        std::vector<Zero> zeroes;
        std::vector<bool> zeroInColumn(matrixSize, false);
        lowerBound = sumCurrentSolution();
       
        createWorkingCopy();

        findRowMinima(rowMinima, zeroInColumn);
        lowerBound += sumMinima(rowMinima);
        subtractMinima(rowMinima, true);
        findColMinima(colMinima, zeroInColumn);
        lowerBound += sumMinima(colMinima);
        subtractMinima(colMinima, false);
        gatherZeroes(rowMinima, colMinima, zeroes);
       
        if (zeroes.empty()) {
            invalidSolution = true;
        }
        else {
            findPenalties(zeroes);
            Zero assignment = assignLargestPenalty(zeroes);

            createLeftNode(solution, assignment);
            int subtourIndex = addToResult(solution, assignment);
            createRightNode(solution, assignment, subtourIndex);
            checkSolutionFound();
           
        }
        free(workingCopy);
    }

    int sumCurrentSolution() {
        int result = 0;
        std::for_each(solution.begin(), solution.end(), [&](std::vector<Selection> v) {
            std::for_each(v.begin(), v.end(), [&](Selection s) {
                result += s.val;
            });
        });
        return result;
    }

    void checkSolutionFound() {
        solutionFound = (solution.size() == 1 && solution[0].size() == matrixSize - 1);
    }

    bool isSolutionFound() {
        return solutionFound;
    }

    std::vector<std::vector<Selection>> getSolution() {
        return solution;
    }

    void createWorkingCopy() {
        workingCopy = (uint16_t*)malloc(sizeof(uint16_t) * matrixSize * matrixSize);
        memcpy(workingCopy, originalMatrix, (sizeof(uint16_t) * matrixSize * matrixSize));
    }

    void findRowMinima(std::vector<Selection>& rowMinima, std::vector<bool>& zeroInColumn) {
        for (int row = 0; row < matrixSize; row++) {
            int rowOffset = row * matrixSize;
            uint16_t currentMinimum = UINT16_MAX;
            Selection currentSelection;
            for (int col = 0; col < matrixSize; col++) {
                if (workingCopy[rowOffset + col] < currentMinimum) {
                    currentMinimum = workingCopy[rowOffset + col];
                    currentSelection = Selection(row, col, currentMinimum);
                    zeroInColumn[col] = true;
                }
            }
            if (currentMinimum != UINT16_MAX) {
                rowMinima.push_back(currentSelection);
            }
        }
        //std::for_each(rowMinima.begin(), rowMinima.end(), [&](Selection s) {
        //    std::cout << " ( " << s.row << ", " << s.col << ") val = " << s.val << " \n";
        //});

    }

    int sumMinima(const std::vector<Selection>& minima) {
        int sum_of_elems = 0;
        std::for_each(minima.begin(), minima.end(), [&](Selection s) {
            sum_of_elems += s.val;
        });
        return sum_of_elems;
    }

    void subtractMinima(std::vector<Selection> minima, const bool row) {
        if (row) {
            std::for_each(minima.begin(), minima.end(), [&](Selection s) {
                if (s.val != 0) {
                    for (int j = 0; j < matrixSize; j++) {
                        if (workingCopy[s.row*matrixSize + j] != UINT16_MAX) {
                            workingCopy[s.row*matrixSize + j] -= s.val;
                        }
                    }
                }
            });
        }
        else {
            std::for_each(minima.begin(), minima.end(), [&](Selection s) {
                if (s.val != 0) {
                    for (int row = 0; row < matrixSize; row++) {
                        if (workingCopy[row*matrixSize + s.col] != UINT16_MAX) {
                            workingCopy[row*matrixSize + s.col] -= s.val;
                        }
                    }
                }
            });
        }
    }

    void findColMinima(std::vector<Selection>& colMinima, const std::vector<bool>& zeroInColumn) {
        for (int col = 0; col < matrixSize; col++) {
            if (!zeroInColumn[col]) {
                uint16_t currentMinimum = UINT16_MAX;
                Selection currentSelection;
                for (int row = 0; row < matrixSize; row++) {
                    if (workingCopy[row*matrixSize + col] < currentMinimum) {
                        currentMinimum = workingCopy[row*matrixSize + col];
                        currentSelection = Selection(row, col, currentMinimum);
                     }
                }
                if (currentMinimum != UINT16_MAX) {
                    colMinima.push_back(currentSelection);
                }
            }
        }
        //std::for_each(colMinima.begin(), colMinima.end(), [&](Selection s) {
        //    std::cout << " ( " << s.row << ", " << s.col << ") val = " << s.val << " \n";
        //});
    }

    void printMatrix() {
        std::cout << "\n Matrix: ";
        for (int i = 0; i < matrixSize; i++) {
            std::cout << "\n";
            for (int j = 0; j < matrixSize; j++) {
                std::cout << workingCopy[i*matrixSize + j] << " ";
            }
        }
        std::cout << "\n";
    }

    void printOriginalMatrix() {
        std::cout << "\nOriginal Matrix: ";
        for (int i = 0; i < matrixSize; i++) {
            std::cout << "\n";
            for (int j = 0; j < matrixSize; j++) {
                std::cout << originalMatrix[i*matrixSize + j] << " ";
            }
        }
        std::cout << "\n";
    }

    void printSolution() {
        std::cout << "\n intermediate :\n";
        std::for_each(solution.begin(), solution.end(), [&](std::vector<Selection> v) {
            std::for_each(v.begin(), v.end(), [&](Selection s) {
                std::cout << "(" << s.row << ", " << s.col << ") value = " << s.val << " ";
            });
            std::cout << "\n";
        });
        std::cout << "\n";
    }
    void gatherZeroes(std::vector<Selection>& rowMinima, std::vector<Selection>& colMinima, std::vector<Zero>& zeroes) {
        for (int i = 0; i < matrixSize; i++) {
            for (int j = 0; j < matrixSize; j++) {
                if (workingCopy[i*matrixSize + j] == 0) {
                    zeroes.push_back(Zero(i, j));
                }
            }
        }
        //std::for_each(zeroes.begin(), zeroes.end(), [&](Zero z) {
        //    std::cout << " zero found : row : " << z.row << " col : " << z.col << "\n";
        //});
    }

    void findPenalties(std::vector<Zero>& zeroes) {
        std::for_each(zeroes.begin(), zeroes.end(), [&](Zero& z) {
            uint16_t horizontalMin = UINT16_MAX;
            for (int col = 0; col < matrixSize; col++) {
                if (col != z.col && workingCopy[z.row*matrixSize + col] < horizontalMin) {
                    horizontalMin = workingCopy[z.row*matrixSize + col];
                }
            }
            z.horPen = horizontalMin;

            uint16_t verticalMin = UINT16_MAX;
            for (int row = 0; row < matrixSize; row++) {
                if (row != z.row && workingCopy[row * matrixSize + z.col] < verticalMin) {
                    verticalMin = workingCopy[row * matrixSize + z.col];
                }
            }
            z.vertPen = verticalMin;
        });

        //std::for_each(zeroes.begin(), zeroes.end(), [&](Zero z) {
        //    std::cout << " (" << z.row << ", " << z.col << ") Horizontal penalty: " << z.horPen << " Vertical penalty: " << z.vertPen << "\n";
        //});
    }

    Zero assignLargestPenalty(std::vector<Zero>& zeroes) {
        Zero assignment = zeroes[0];
        std::for_each(zeroes.begin(), zeroes.end(), [&](Zero z) {
            if ((z.vertPen + z.horPen) > (assignment.vertPen + assignment.horPen)) {
                assignment = z;
            }
        });
        return assignment;
    }

    int addToResult(std::vector<std::vector<Selection>>& solution, Zero assignment) {
        int startOfSubtourIndex = -1;
        int endOfSubtourIndex = -1;
        int index = 0;
        int result = 0;

        //dealing with subtours....
        std::for_each(solution.begin(), solution.end(), [&](std::vector<Selection>& subtour) {
            if (subtour[0].row == assignment.col) {
                startOfSubtourIndex = index;
            }
            else if (subtour.back().col == assignment.row) {
                endOfSubtourIndex = index;
            }
            index++;
        });
        if (startOfSubtourIndex > -1 && endOfSubtourIndex > -1) {
            //this case is true if our solution connects 2 subtours, by being at the start of one and at the end of the other.
            //we combine them by concatenating the subtour that it is the start of after the subtour that it is the end of.
            solution[endOfSubtourIndex].push_back(Selection(assignment.row, assignment.col, originalMatrix[assignment.row* matrixSize + assignment.col]));
            solution[endOfSubtourIndex].insert(solution[endOfSubtourIndex].end(), solution[startOfSubtourIndex].begin(), solution[startOfSubtourIndex].end());
            solution.erase(solution.begin() + startOfSubtourIndex);
            result = endOfSubtourIndex;
            if (startOfSubtourIndex < endOfSubtourIndex) {
                result--;
            }

        }
        else if (startOfSubtourIndex > -1) {
            solution[startOfSubtourIndex].insert(solution[startOfSubtourIndex].begin(), Selection(assignment.row, assignment.col, originalMatrix[assignment.row * matrixSize + assignment.col]));
            //subtour elimination, we can no longer go FROM(column) subtour.back().row TO(row) assignment.col.
            result = startOfSubtourIndex;
        }
        else if (endOfSubtourIndex > -1) {
            solution[endOfSubtourIndex].push_back(Selection(assignment.row, assignment.col, originalMatrix[assignment.row* matrixSize + assignment.col]));
            //subtour elimination, we can no longer go FROM(column) assignment.row TO(row) subtour[0].col.
            result = endOfSubtourIndex;
        }
        else {
            solution.push_back(std::vector<Selection>(1, Selection(assignment.row, assignment.col, originalMatrix[assignment.row* matrixSize + assignment.col])));
            result = solution.size() - 1;
        }
        return result;
    }

    void createLeftNode(std::vector<std::vector<Selection>>& solution, Zero assignment) {
        uint16_t* leftSideMatrix = (uint16_t*)malloc(sizeof(uint16_t)*matrixSize*matrixSize);
        for (int row = 0; row < matrixSize; row++) {
            for (int col = 0; col < matrixSize; col++) {
                if (row == assignment.row && col == assignment.col) {
                    leftSideMatrix[row*matrixSize + col] = UINT16_MAX;
                }
                else {
                    leftSideMatrix[row*matrixSize + col] = originalMatrix[row*matrixSize + col];
                }
            }
        }
        left = new Node(matrixSize, leftSideMatrix, solution, lowerBound + assignment.horPen + assignment.vertPen);
    }

    void createRightNode(std::vector<std::vector<Selection>>& solution, Zero assignment, int subtourIndex) {
        uint16_t* rightSideMatrix = (uint16_t*)malloc(sizeof(uint16_t)*matrixSize*matrixSize);
        for (int row = 0; row < matrixSize; row++) {
            for (int col = 0; col < matrixSize; col++) {
                if (assignment.row == row || assignment.col == col) {
                    rightSideMatrix[row*matrixSize + col] = UINT16_MAX;
                }
                else {
                    rightSideMatrix[row*matrixSize + col] = originalMatrix[row*matrixSize + col];
                }
            }
        }
        rightSideMatrix[solution[subtourIndex].back().col * matrixSize + solution[subtourIndex][0].row] = UINT16_MAX;
        right = new Node(matrixSize, rightSideMatrix, solution);
    }

    Node getRight() {
        return *right;
    }

    Node getLeft() {
        return *left;
    }

};

Node::Node(int amountOfCities, uint16_t* matrix, std::vector<std::vector<Selection>> solution) : solution(solution) {
    matrixSize = amountOfCities;
    originalMatrix = matrix;
}

Node::Node(int amountOfCities,uint16_t* matrix, std::vector<std::vector<Selection>> solution, uint16_t lb) : solution(solution){
    matrixSize = amountOfCities;
    originalMatrix = matrix;
    lowerBound = lb;
}



class Tree {
    Node *start;
    int lowerBound;
    int upperBound;
    std::vector<Node> nodes;
};

void getMatrixFromFile(std::string filename, int matrixSize, uint16_t* originalMatrix) {
    std::vector<int> xValues(matrixSize);
    std::vector<int> yValues(matrixSize);
    int n, m, i = 0;
    std::ifstream read(filename);
    while (i < matrixSize && read >> n >> m) {
        xValues[i] = n;
        yValues[i] = m;
        i++;
    }
   
    for (int i = 0; i < matrixSize; i++) {
        for (int j = 0; j < matrixSize; j++) {
            if (i == j) {
                originalMatrix[i*matrixSize + j] = UINT16_MAX;
            }
            else {
                originalMatrix[i*matrixSize + j] = (uint16_t)sqrt(pow(xValues[i]-xValues[j], 2) + pow(yValues[i] - yValues[j],2));
            }
        }
    }
}


void pruneNodeList(std::vector<Node>& nodeList, int upperBound) {
    int i = 0;
    while (i < nodeList.size()) {
        if (nodeList[i].lowerBound >= upperBound) {
            nodeList.erase(nodeList.begin() + i);
        }
        else {
            ++i;
        }
    }
}

int calculateSolution(uint16_t* originalMatrix, std::vector<std::vector<Selection>>& solution, int matrixSize) {
    int sum = 0;
    std::for_each(solution.begin(), solution.end(), [&](std::vector<Selection> v) {
        std::for_each(v.begin(), v.end(), [&](Selection s) {
            sum += s.val;
        });
    });
    sum += originalMatrix[solution[0].back().col * matrixSize + solution[0][0].row];
    return sum;
}

int findNextEvaluation(std::vector<Node>& nodeList) {
    Node eval = nodeList[0];
    int index = 0;
    int currentIndex = 0;
    std::for_each(nodeList.begin(), nodeList.end(), [&](Node n) {
        if (n.lowerBound < eval.lowerBound) {
            currentIndex = index;
        }
        index++;
    });
    return currentIndex;
}

int main() {
    int amountOfCities = 15;
    uint16_t* originalMatrix = (uint16_t*)malloc(sizeof(uint16_t)*amountOfCities*amountOfCities);
    getMatrixFromFile("cities.txt", amountOfCities, originalMatrix);
    std::vector<std::vector<Selection>> solution;
    Node node = Node(amountOfCities, originalMatrix, solution);
    std::vector<Node> nodeList;

    std::vector<std::vector<Selection>> currentSolution;
    std::vector<std::vector<Selection>> tempSolution;
    int upperBound = UINT16_MAX;
    int count = 0;
    clock_t t1, t2;
    t1 = clock();
    while (!node.isSolutionFound()) {
        node.evaluate();
        count++;
        if (!node.isSolutionFound()) {
            nodeList.push_back(node.getLeft());
            node = node.getRight();
        }
    }
    currentSolution = node.getSolution();
    upperBound = calculateSolution(originalMatrix, currentSolution, amountOfCities);
    pruneNodeList(nodeList, upperBound);
    while (!nodeList.empty()) {
        int nextNode = findNextEvaluation(nodeList);
        node = nodeList[nextNode];
        nodeList.erase(nodeList.begin() + nextNode);
        while (!node.isSolutionFound() && !node.invalidSolution) {
            node.evaluate();
            count++;
            if (!node.isSolutionFound() && !node.invalidSolution && node.lowerBound < upperBound) {
                if (node.getLeft().lowerBound < upperBound) {
                    nodeList.push_back(node.getLeft());
                }
                node = node.getRight();
            }
            else if (node.lowerBound >= upperBound) {
                break;
            } else if(!node.invalidSolution){
                tempSolution = node.getSolution();
                int value = calculateSolution(originalMatrix, tempSolution, amountOfCities);
                if (value < upperBound) {
                    currentSolution = tempSolution;
                    upperBound = value;
                    pruneNodeList(nodeList, upperBound);
                }
            }
        }

    }
    t2 = clock();
    float diff((float)t2 - (float)t1);
    std::cout << "time spent: " << diff << "\n";
    std::cout << "count = " << count << "\n";
    std::cout << "\n Solution :\n";
    std::for_each(currentSolution.begin(), currentSolution.end(), [&](std::vector<Selection> v) {
        std::for_each(v.begin(), v.end(), [&](Selection s) {
            std::cout << "(" << s.row << ", " << s.col << ") ";
        });
        std::cout << "\n";
    });
    std::cout << "\nResult = " << upperBound;
}