<file_sep>#ifndef KDTREE_HPP_
#define KDTREE_HPP_

#include <algorithm>
#include <cfloat>
#include <vector>
#include "Point.hpp"

template <size_t D>
class KDTree {
 private:
  class KDTreeNode {
   public:
    /* Constructors */
    KDTreeNode() : KDTreeNode{0, 0} {}
    KDTreeNode(size_t axis) : KDTreeNode{axis, 0} {}
    KDTreeNode(size_t axis, double value)
        : KDTreeNode{{}, axis, value, nullptr, nullptr} {}
    KDTreeNode(size_t axis, Point<D> point)
        : KDTreeNode{point, axis, point[axis], nullptr, nullptr} {}
    KDTreeNode(const Point<D>& point, size_t axis, double value,
               KDTreeNode* left, KDTreeNode* right)
        : point_{point},
          axis_{axis},
          value_{value},
          left_{left},
          right_{right} {}
    /**
     * Deep destructor
     */
    ~KDTreeNode() {
      if (left_ != nullptr) delete left_;
      if (right_ != nullptr) delete right_;
    }
    bool isLeaf() const { return left_ == nullptr && right_ == nullptr; }
    Point<D> point_;
    size_t axis_;
    double value_;
    KDTreeNode* left_;
    KDTreeNode* right_;
  };

 public:
  KDTreeNode* root;

  KDTree(std::vector<Point<D>> points) : root{createTree(points, 0)} {}
  ~KDTree() {
    if (root != nullptr) delete root;
  }

  void deleteNode(KDTreeNode* node) {
    if (node != nullptr) delete node;
  }

  KDTreeNode* createTree(std::vector<Point<D>> points, size_t depth) {
    size_t axis = depth % D;
    std::sort(points.begin(), points.end(),
              [axis](const Point<D>& a, const Point<D>& b) {
                return a[axis] < b[axis];
              });
    if (points.size() == 1) {
      return new KDTreeNode{axis, points[0]};
    } else if (points.size() < 1) {
      return nullptr;
    }
    size_t middleIndex = points.size() / 2;  // For indexing points
    int midInc = middleIndex;                // For adding to iterator
    KDTreeNode* node = new KDTreeNode{axis, points[middleIndex]};
    node->left_ =
        createTree({points.begin(), points.begin() + midInc}, depth + 1);
    node->right_ =
        createTree({points.begin() + midInc, points.end()}, depth + 1);
    return node;
  }

  Point<D> nns(const Point<D>& point) {
    double infinity = DBL_MAX;
    Point<D> nearestNeighbor{};
    nnsRecurse(point, root, infinity, nearestNeighbor);
    return nearestNeighbor;
  }

  void nnsRecurse(const Point<D>& point, KDTreeNode* node, double& radius,
                  Point<D>& nearestNeighbor) {
    if (node == nullptr) return;
    if (node->isLeaf()) {
      double pointsDistance = distance(point, node->point_);
      if (radius > pointsDistance) {
        radius = pointsDistance;
        nearestNeighbor = node->point_;
      }
      return;
    }
    if (point[node->axis_] < node->value_) {
      nnsRecurse(point, node->left_, radius, nearestNeighbor);
      if (radius == DBL_MAX || point[node->axis_] + radius >= node->value_) {
        nnsRecurse(point, node->right_, radius, nearestNeighbor);
      }
    } else {
      nnsRecurse(point, node->right_, radius, nearestNeighbor);
      if (radius == DBL_MAX || point[node->axis_] - radius < node->value_) {
        nnsRecurse(point, node->left_, radius, nearestNeighbor);
      }
    }
  }
};

#endif  // KDTREE_HPP_
<file_sep>CXX ?= g++
CXXFLAGS := -g -O3 -std=c++14 -Wall -Wextra -pedantic
LDFLAGS :=

all: main.o NearestNeighbor.o
	$(CXX) $(LDFLAGS) main.o NearestNeighbor.o -o main.exe

main.o: main.cpp KDTree.hpp Point.hpp
	$(CXX) $(CXXFLAGS) -c main.cpp

NearestNeighbor.o: NearestNeighbor.cpp NearestNeighbor.hpp KDTree.hpp Point.hpp
	$(CXX) $(CXXFLAGS) -c NearestNeighbor.cpp

clean:
	-rm *.o main.exe
<file_sep>#include "NearestNeighbor.hpp"
#include <chrono>
#include <fstream>
#include <iostream>
#include <stdexcept>

#define DIMENSION 3
#define NUM_POINTS 3e6

void run_tests() {
  std::cout << "# of points: " << NUM_POINTS << std::endl;
  auto cloud = generate_point_cloud<DIMENSION>(NUM_POINTS, {-100.0, 100.0});
  auto point = random_point<DIMENSION>({-100.0, 100.0});
  std::cout << "---------- Brute Force Method ---------- " << std::endl;
  auto v = makeDecorator(brute_force<DIMENSION>)(point, cloud, 1);
  std::cout << "Closest points to " << point << std::endl;
  for (auto e : v) {
    std::cout << e << std::endl;
    std::cout << "Distance: " << distance(point, e) << std::endl;
  }
  std::cout << "---------- KD-Tree Method ---------- " << std::endl;
  std::cout << "> Constructing tree..." << std::endl;
  auto kdTree = makeDecorator(initialize_kd_tree<DIMENSION>)(cloud);
  std::cout << "> Running NNS algorithm..." << std::endl;
  auto nearestNeighbor = makeDecorator(kd_tree_nns<DIMENSION>)(point, kdTree);
  std::cout << "Closest point to " << point << std::endl;
  std::cout << nearestNeighbor << std::endl;
  std::cout << "Distance: " << distance(point, nearestNeighbor) << std::endl;

  std::ofstream file{"results.csv"};
  if (!file.is_open()) throw std::runtime_error{"Unable to open file"};
  file << "Points,Construction Time (ns),NNS Time (ns)" << std::endl;
  file << std::fixed;
  for (int i = 0; i < 8; i++) {
    std::cout << "10^" << i << " points" << std::endl;
    cloud = generate_point_cloud<DIMENSION>(pow(10, i), {-100.0, 100.0});
    point = random_point<DIMENSION>({-100.0, 100.0});
    std::chrono::time_point<std::chrono::system_clock> start, end;
    start = std::chrono::system_clock::now();
    initialize_kd_tree<DIMENSION>(cloud);
    end = std::chrono::system_clock::now();
    std::chrono::nanoseconds kd_time = end - start;

    start = std::chrono::system_clock::now();
    kd_tree_nns<DIMENSION>(point, kdTree);
    end = std::chrono::system_clock::now();
    std::chrono::nanoseconds nns_time = end - start;
    file << static_cast<int>(pow(10, i)) << "," << kd_time.count() << ","
         << nns_time.count() << std::endl;
  }
}
<file_sep># cs375-final-project
CS 375 Final Project - Nearest Neighbor
<file_sep>#include <vector>
#include "KDTree.hpp"
#include "NearestNeighbor.hpp"
#include "Point.hpp"

int main() {
  run_tests();
  return 0;
}
<file_sep>cmake_minimum_required(VERSION 3.2 FATAL_ERROR)
project(NearestNeighbor)
set(CMAKE_BUILD_TYPE Release)

## C++11 support
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

## Increase warning levels
if(CMAKE_COMPILER_IS_GNUCXX)
    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pedantic -Wall -Wextra -Wcast-align -Wcast-qual -Wctor-dtor-privacy -Wdisabled-optimization -Wformat=2 -Winit-self -Wlogical-op -Wmissing-declarations -Wmissing-include-dirs -Wnoexcept -Wold-style-cast -Woverloaded-virtual -Wredundant-decls -Wshadow -Wsign-conversion -Wsign-promo -Wstrict-null-sentinel -Wstrict-overflow=5 -Wswitch-default -Wundef -Werror")
    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g")
endif()

add_executable(main.exe main.cpp NearestNeighbor.cpp)
target_compile_features(main.exe PRIVATE
    cxx_aggregate_default_initializers
    cxx_lambdas
    cxx_nullptr
    cxx_auto_type
    cxx_generalized_initializers
    cxx_range_for
    cxx_uniform_initialization
    cxx_delegating_constructors)
<file_sep>#ifndef NEARESTNEIGHBOR_HPP_
#define NEARESTNEIGHBOR_HPP_

#include <algorithm>
#include <chrono>
#include <functional>
#include <random>
#include <utility>
#include <vector>
#include "KDTree.hpp"

static std::random_device rand_dev;

template <typename>
struct Decorator;

void run_tests();

/*
 * For creating a single point.
 * Not used in generate_point_cloud because of efficiency
 */
template <size_t D>
Point<D> random_point(const std::pair<double, double>& range) {
  std::default_random_engine e1(rand_dev());
  std::uniform_real_distribution<double> uniform_dist{range.first,
                                                      range.second};
  std::array<double, D> components;
  for (size_t j = 0; j < D; j++) {
    components[j] = uniform_dist(e1);
  }
  return Point<D>{components};
}

template <size_t D>
std::vector<Point<D>> generate_point_cloud(
    size_t number, const std::pair<double, double>& range) {
  std::default_random_engine e1(rand_dev());
  std::uniform_real_distribution<double> uniform_dist{range.first,
                                                      range.second};
  std::vector<Point<D>> ret;
  for (size_t i = 0; i < number; i++) {
    std::array<double, D> components;
    for (size_t j = 0; j < D; j++) {
      components[j] = uniform_dist(e1);
    }
    ret.emplace_back(components);
  }
  return ret;
}

template <size_t D>
std::vector<Point<D>> brute_force(const Point<D>& point,
                                  std::vector<Point<D>> neighbors, size_t num) {
  std::vector<Point<D>> nearest;
  for (size_t i = 0; i < num; i++) {
    auto neighbor =
        std::min_element(neighbors.begin(), neighbors.end(),
                         [point](const Point<D>& lhs, const Point<D>& rhs) {
                           return distance(point, lhs) < distance(point, rhs);
                         });
    nearest.push_back(*neighbor);
    neighbors.erase(neighbor);
  }
  return nearest;
}

template <size_t D>
Point<D> kd_tree_nns(const Point<D>& point, KDTree<D>& kdTree) {
  return kdTree.nns(point);
}

template <size_t D>
KDTree<D> initialize_kd_tree(const std::vector<Point<D>>& neighbors) {
  return KDTree<D>{neighbors};
}

template <size_t D>
Point<D> kd_nearest_neighbor(const Point<D>& point, const KDTree<D>& kdTree) {
  Point<D> nearestNeighbor = kd_tree_nns(point, kdTree);
  return nearestNeighbor;
}

/* http://stackoverflow.com/questions/30679445/python-like-c-decorators */
template <typename R, typename... Args>
struct Decorator<R(Args...)> {
  Decorator(std::function<R(Args...)> f) : f_(f) {}

  R operator()(Args... args) {
    std::chrono::time_point<std::chrono::system_clock> start, end;
    start = std::chrono::system_clock::now();
    R ret = f_(args...);
    end = std::chrono::system_clock::now();
    std::chrono::nanoseconds elapsed_nanos =
        std::chrono::duration_cast<std::chrono::nanoseconds>(end - start);
    std::cout << "Execution took " << elapsed_nanos.count() << "ns"
              << std::endl;
    return ret;
  }
  std::function<R(Args...)> f_;
};

template <typename R, typename... Args>
Decorator<R(Args...)> makeDecorator(R (*f)(Args...)) {
  return Decorator<R(Args...)>(std::function<R(Args...)>(f));
}
#endif  // NEARESTNEIGHBOR_HPP_
