From 24ee655e6823bb2d83171638e2c53166ea71f6bd Mon Sep 17 00:00:00 2001 From: Garry Boyer Date: Wed, 14 Mar 2007 17:18:28 +0000 Subject: [PATCH] hi --- fastlib/fastlib/build.py | 10 + fastlib/file/build.py | 7 + fastlib/file/serialize.cc | 79 ++++++++ fastlib/file/serialize.h | 237 +++++++++++++++++++++++ fastlib/par/build.py | 12 ++ fastlib/par/grain.h | 213 ++++++++++++++++++++ fastlib/par/mpigrain.h | 398 ++++++++++++++++++++++++++++++++++++++ fastlib/par/task.h | 30 +++ fastlib/par/thread.cc | 5 + fastlib/par/thread.h | 203 +++++++++++++++++++ fastlib/tree/bounds.h | 353 +++++++++++++++++++++++++++++++++ fastlib/tree/build.py | 10 + fastlib/tree/kdtree.h | 253 ++++++++++++++++++++++++ fastlib/tree/spacetree.h | 279 ++++++++++++++++++++++++++ fastlib/tree/statistic.h | 40 ++++ 15 files changed, 2129 insertions(+) create mode 100644 fastlib/file/serialize.cc create mode 100644 fastlib/file/serialize.h create mode 100644 fastlib/par/build.py create mode 100644 fastlib/par/grain.h create mode 100644 fastlib/par/mpigrain.h create mode 100644 fastlib/par/task.h create mode 100644 fastlib/par/thread.cc create mode 100644 fastlib/par/thread.h create mode 100644 fastlib/tree/bounds.h create mode 100644 fastlib/tree/build.py create mode 100644 fastlib/tree/kdtree.h create mode 100644 fastlib/tree/spacetree.h create mode 100644 fastlib/tree/statistic.h diff --git a/fastlib/fastlib/build.py b/fastlib/fastlib/build.py index 8bcbd48cde..510eb32e6f 100644 --- a/fastlib/fastlib/build.py +++ b/fastlib/fastlib/build.py @@ -7,3 +7,13 @@ librule( #, "tree:tree", "par:par" ] ) + +librule( + name = "fastlib_int", + headers = ["fastlib.h"], + deplibs = ["la:la", "base:base", + "fx:fx", "file:file", "col:col", + "data:data", "math:math" + "tree:tree", "par:par" + ] + ) diff --git a/fastlib/file/build.py b/fastlib/file/build.py index ca7bf245d2..4aae2ce1b2 100644 --- a/fastlib/file/build.py +++ b/fastlib/file/build.py @@ -5,6 +5,13 @@ librule( deplibs = ["base:base", "col:col"], ) +librule( # internal rule + name = "file_int", + sources = ["serialize.cc"], + headers = ["serialize.h"], + deplibs = [":file"], + ) + binrule( name = "textfile_test", sources = ["textfile_test.cc"], diff --git a/fastlib/file/serialize.cc b/fastlib/file/serialize.cc new file mode 100644 index 0000000000..66b149e5dc --- /dev/null +++ b/fastlib/file/serialize.cc @@ -0,0 +1,79 @@ +/** + * @file serialize.cc + * + * Definitions for serialization methods. + */ + +#include "serialize.h" + +#include +#include +#include +#include + +file::magic_t file::CreateMagic(const char *str) { + magic_t magic = 0x31415926; + + for (; *str != '\0'; str++) { + magic = (magic << 9) + (magic >> 11) + *str; + } + + return magic; +} + +success_t NativeArraySerializer::WriteFile(const char *fname) { + FILE *f = fopen(fname, "wb"); + + if (unlikely(!f)) { + return SUCCESS_FAIL; + } + + if (unlikely(index_t(fwrite(ptr(), 1, size(), f)) != index_t(size()))) { + (void) fclose(f); + return SUCCESS_FAIL; + } + + return SUCCESS_FROM_INT(fclose(f)); +} + +success_t NativeFileDeserializer::Init(const char *fname) { + struct stat info; + FILE *f; + size_t file_size; + char *data; + + DEBUG_MSG(1.0, "Opening a native file deserializer."); + + if (unlikely(stat(fname, &info) < 0)) { + NONFATAL("File [%s] does not exist.", fname); + data = NULL; + } else { + f = fopen(fname, "rb"); + + if (unlikely(!f)) { + data = NULL; + NONFATAL("Cannot read [%s], but it exists.", fname); + } else { + file_size = size_t(info.st_size); + data = mem::Alloc(file_size); + + if (unlikely(fread(data, 1, file_size, f) != file_size)) { + mem::Free(data); + data = NULL; + NONFATAL("Error reading contents of [%s].", fname); + } + + (void) fclose(f); + } + } + + if (data == NULL) { + // We have to initialize this into a valid state, or else the program + // will segfault. + file_size = 0; + } + + NativeArrayDeserializer::Init(data, file_size); + + return (data != NULL) ? SUCCESS_PASS : SUCCESS_FAIL; +} diff --git a/fastlib/file/serialize.h b/fastlib/file/serialize.h new file mode 100644 index 0000000000..8f22fb4748 --- /dev/null +++ b/fastlib/file/serialize.h @@ -0,0 +1,237 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @file serialize.h + * + * Tools for serialization. + * + * See NativeArraySerializer and NativeArrayDeserializer for an example of + * the Serializer and Deserializer "concepts". + * + * WARNING: This has been only casually tested, not thoroughly tested. + * Use at your own risk. + */ + +#ifndef FILE_SERIALIZE_H +#define FILE_SERIALIZE_H + +#include "base/common.h" +#include "base/ccmem.h" +#include "col/arraylist.h" + +#include + +namespace file { + /** + * Magic number type. + */ + typedef uint32 magic_t; + + /** + * Automatically creates a magic number by hashing the + * string. + */ + magic_t CreateMagic(const char *str); +}; + +/** + * Automatically creates a magic number for a class. + * + * This is especially useful if you have a templated class. + * You can add other magic numbers to it if you want the number to be + * more specific, like MAGIC_NUMBER(a) + MAGIC_NUMBER(b). + * + * @param classname the name of the class + */ +#define MAGIC_NUMBER(classname) \ + (file::CreateMagic(typeid(classname).name())) + +/** + * Use this to serialize a stream of objects in an array. + * + * Useful before sending canonical versions over MPI, or to a file. + * + * The computer reading this should use the same CPU type. + * In particular, the type sizes and byte order must remain the same. + * + * Only primitives or structs of primitives should be stored; in + * particular, no pointers. + * + * TODO: More complicated serializers (like text serializers) later + * might have a slight problem with serializing structs. Come back + * later for info. Maybe we'll decide against text serializers. + */ +class NativeArraySerializer { + FORBID_COPY(NativeArraySerializer); + + private: + ArrayList data_; + + public: + NativeArraySerializer() {} + ~NativeArraySerializer() {} + + void Init() { + data_.Init(); + } + + void PutMagic(file::magic_t magic_num) { + Put(magic_num); + } + + /** + * Appends a struct or primitive to the end. + * + * This will store the object bit-by-bit, so don't store any + * pointers. + */ + template + void Put(const T& val) { + mem::CopyBytes(data_.AddBack(sizeof(T)), &val, sizeof(T)); + } + + /** + * Appends an array of structs or primtives to the end. + * + * This will store the object bit-by-bit, so don't store any + * pointers. + */ + template + void Put(const T* array, index_t count) { + size_t bytes = count * sizeof(T); + mem::CopyBytes(data_.AddBack(bytes), array, bytes); + } + + /** + * Returns the data serialized, clearing the internal + * state of this serializer. + * + * The pointer must eventually be freed with mem::Free. + * + * @return a pointer that can later be passed to deserializer + */ + char *ReleasePointer() { + data_.Trim(); + return data_.ReleasePointer(); + } + + const char *ptr() const { + return data_.begin(); + } + + index_t size() const { + return data_.size(); + } + + /** + * Dumps to a file. + */ + success_t WriteFile(const char *fname); +}; + +/** + * Use this to deserialize a stream created by NativeArraySerializer. + * + * This computer should use the same architecture as was used to seralize. + * Specifically, the type sizes and byte order must be identical. + */ +class NativeArrayDeserializer { + private: + const char *ptr_; + index_t pos_; + index_t size_; + + public: + NativeArrayDeserializer() { + DEBUG_POISON_PTR(ptr_); + } + ~NativeArrayDeserializer() {} + + /** + * Initializes given a pointer. + * + * The pointer will not be freed. + */ + void Init(const char *ptr_in, index_t size_in) { + ptr_ = ptr_in; + size_ = size_in; + pos_ = 0; + } + + success_t CheckMagic(file::magic_t magic_num) { + file::magic_t x; + Get(&x); + return likely(x == magic_num) ? SUCCESS_PASS : SUCCESS_FAIL; + } + + void AssertMagic(file::magic_t magic_num) { + file::magic_t x; + Get(&x); + assert(x == magic_num); + } + + /** + * Retrieves the next struct or primtiive. + */ + template + void Get(T *dest) { + DEBUG_ASSERT(size_t(pos_ + sizeof(T)) <= size_t(size_)); + mem::CopyBytes(dest, ptr_ + pos_, sizeof(T)); + pos_ += sizeof(T); + } + + /** + * Retrieves the next array of structs or primitives. + */ + template + void Get(T *dest, index_t count) { + size_t bytes = count * sizeof(T); + DEBUG_ASSERT(index_t(pos_ + bytes) <= index_t(size_)); + mem::CopyBytes(dest, ptr_ + pos_, bytes); + pos_ += bytes; + } + + /** + * Returns whether finished. + * + * You are not recommended to use this; store counts instead. + * + * @return true iff there is no data remaining + */ + bool done() const { + return pos_ == size_; + } + + size_t size() const { + return size_; + } + + size_t pos() const { + return pos_; + } + + const char *ptr() const { + return ptr_; + } +}; + +/** + * Deserializer from a file. + * + * NOTE: This currently uses a NativeArrayDeserializer, but in the future + * may use buffering. + */ +class NativeFileDeserializer : public NativeArrayDeserializer { + public: + NativeFileDeserializer() {} + ~NativeFileDeserializer() { + mem::Free(const_cast(ptr())); + } + + /** + * Initializes by reading the entire file into memory. + */ + success_t Init(const char *fname); +}; + +#endif diff --git a/fastlib/par/build.py b/fastlib/par/build.py new file mode 100644 index 0000000000..e3318b546f --- /dev/null +++ b/fastlib/par/build.py @@ -0,0 +1,12 @@ + +librule( + sources = ["thread.cc"], + headers = ["thread.h", "task.h", "grain.h"], + deplibs = ["base:base", "col:col"]) + +librule( + name = "mpi", + sources = [], + headers = ["mpigrain.h"], + deplibs = [":par"]) + diff --git a/fastlib/par/grain.h b/fastlib/par/grain.h new file mode 100644 index 0000000000..0271482229 --- /dev/null +++ b/fastlib/par/grain.h @@ -0,0 +1,213 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @file grain.h + * + * Tool for creating very simple parallel programs by enqueuing small work + * items in a priority queue and executing them greedily. + */ + +#ifndef PAR_GRAIN_H +#define PAR_GRAIN_H + +#include "col/heap.h" + +#include "thread.h" + +/** + * Simple difficulty-based work queue for easy parallelization. + * + * To use this, simply divide up your work into grains (as shown below) that + * are probably a bit smaller than what one thread should be able to handle. + * Enqueue each grain into the GrainQueue, associating it with a + * "difficulty" measure which should estimate the relative amount of time + * for each grain. You can then use this to automatically run a number of + * threads. + * + * To allow grains to be sent over multiple machines, put only basic data + * into your grains, and associate the grain queue with a context, which + * is a pointer to something. This context will be passed to every grain + * when it is run. + * + * TODO: This has a bug that, if any of the threads find an empty queue, + * that thread will die. This would prohibit you from recursively building + * a kd-tree or such, as after the root node is de-queued all the other + * threads would die because they think there is no work to do. + * + * This class is thread safe, so it is perfectly fine for grains to put + * more work on the grain queue. + * + * TODO: Update documentation to reflect the fact that ThreadedGrainRunner + * is a separate class now. + * + * @code + * struct SolverGrain { + * Solver *solver; + * int a; + * int b; + * ~SolverGrain() {} + * SolverGrain(Solver *solver_in, int a_in, int b_in) { + * solver = solver_in; + * a = a_in; + * b = b_in; + * } + * void Run() { + * solver->Solve(a, b); + * } + * }; + * + * class Solver { + * void SolveRange(int a_min, int a_max, int b_min, int b_max) { + * GrainQueue<SolverGrain> queue; + * for (int a = 0; a < 10; a++) { + * for (int b = 0; b < 10; b++) { + * queue->Put(a * b, new SolverGrain(this, a, b)); + * } + * } + * } + * void Solve(int a, int b) {....} + * } + * @endcode + */ +template +class GrainQueue { + public: + typedef TGrain Grain; + + private: + MinHeap queue_; + Mutex mutex_; + + public: + GrainQueue() {} + ~GrainQueue() {} + + /** + * Initializes. + */ + void Init() { + queue_.Init(); + } + + /** + * Puts a grain into the queue to be dispatched. + * + * @param difficulty relative problem difficulty + */ + void Put(double difficulty, Grain *grain) { + mutex_.Lock(); + queue_.Put(-difficulty, grain); + mutex_.Unlock(); + } + + /** + * Pops the most desirable grain to work on. + * + * You might not have to call this yourself. + */ + Grain *Pop() { + mutex_.Lock(); + Grain *result = likely(queue_.size() != 0) ? queue_.Pop() : NULL; + mutex_.Unlock(); + return result; + } + + /** + * Gets the size of this queue. + */ + index_t size() const { + return queue_.size(); + } +}; + +template +class ThreadedGrainRunner { + FORBID_COPY(ThreadedGrainRunner); + public: + typedef TGrain Grain; + typedef TContext Context; + + private: + struct ThreadTask : public Task { + ThreadedGrainRunner *runner_; + + ThreadTask(ThreadedGrainRunner *runner_in) { + runner_ = runner_in; + } + + void Run() { + while (runner_->RunOneGrain()) {} + delete this; + } + }; + + private: + GrainQueue *queue_; + Context context_; + + public: + ThreadedGrainRunner() {} + ~ThreadedGrainRunner() {} + + void Init(GrainQueue *queue_in, Context context_in) { + queue_ = queue_in; + context_ = context_in; + } + + /** + * Pops and runs one task. + * + * Use this if you, for some reason, decided that running separate threads + * was a bad idea and you really just want to run grains yourself. + * + * @return true whether a task was run, false if no more tasks left + */ + bool RunOneGrain() { + Grain *grain = queue_->Pop(); + if (unlikely(!grain)) { + return false; + } else { + grain->Run(context_); + delete grain; + return true; + } + } + + /** + * Spawns a single running thread. + * + * You must WaitStop or Detach this thread, and eventually, free + * the returned object. + * + * (TODO: In the future you might have to delete its task() too). + * + * @return a newly created thread + */ + Thread *SpawnThread() { + ThreadTask *task = new ThreadTask(this); + Thread *thread = new Thread(); + thread->Init(task); + thread->Start(); + return thread; + } + + /** + * Creates the specified number of threads, and uses those to execute + * all grains of work. + */ + void RunThreads(int num_threads) { + ArrayList threads; + + threads.Init(num_threads); + + for (int i = 0; i < num_threads; i++) { + threads[i] = SpawnThread(); + } + for (int i = 0; i < num_threads; i++) { + threads[i]->WaitStop(); + delete threads[i]; + } + } +}; + +#endif diff --git a/fastlib/par/mpigrain.h b/fastlib/par/mpigrain.h new file mode 100644 index 0000000000..81bd87a6b8 --- /dev/null +++ b/fastlib/par/mpigrain.h @@ -0,0 +1,398 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @file mpigrain.h + * + * Tool for creating simple MPI parallel programs. + * You must build with --compiler=mpi. + * + * WARNING! This currently DOES NOT have ANY WAY for you to send + * results back! + * + * TODO: Broken + */ + +#ifndef PAR_MPIGRAIN_H +#define PAR_MPIGRAIN_H + +#include "xrun/xrun.h" + +#include "grain.h" + +#include + +/** + * Grain-based parallelism over MPI. + * + * Grains must be bit-copiable. In the future serialization might be a better + * idea. + * + * Uses single master that maintains priority queue, and handles work out + * to slaves. + * + * TODO: This does NOT have any way for you to send results back to the + * mater computer. + */ +template +class MPIGrainRunner { + FORBID_COPY(MPIGrainRunner); + friend class MPIDispatcher; + + public: + typedef TGrain Grain; + typedef TContext Context; + + class Dispatcher { + FORBID_COPY(Dispatcher); + friend class MPIGrainRunner; + + private: + class MPIMasterTask : public Task { + private: + struct MPIGrainRunner *runner_; + + public: + MPIMasterTask(MPIGrainRunner *runner_in) { + runner_ = runner_in; + } + + void Run() { + int rank; + int message; + int n_slaves_alive = 0; + int n_slaves_busy = 0; + + DEBUG_MSG(1.0, "DISPATCH: Firing up the cannons."); + DEBUG_MSG(1.0, "DISPATCH: We will accomplish %u tasks.", + unsigned(runner_->dispatcher_->queue_->size())); + + for (;;) { + runner_->RecvInt_(&rank, &message); + + runner_->mutex_.Lock(); + Grain* grain = runner_->dispatcher_->queue_->Pop(); + runner_->mutex_.Unlock(); + + if (message == BIRTH) { + DEBUG_MSG(1.0, "DISPATCH: Received birth message."); + n_slaves_alive++; + n_slaves_busy++; + } else { + DEBUG_ASSERT(message == GIVE_ME_WORK); + DEBUG_MSG(1.0, "DISPATCH: Somebody wants more work."); + } + + if (grain) { + char buf[sizeof(Grain) + 1]; + + buf[0] = 1; // data is available + mem::CopyBytes(buf + 1, grain, sizeof(*grain)); + delete grain; + + DEBUG_MSG(1.0, "DISPATCH: Sending work on over."); + MPI_Send(buf, sizeof(buf), MPI_CHAR, rank, + runner_->tag_, MPI_COMM_WORLD); + } else if (n_slaves_alive == runner_->n_slaves_) { + // we have to make sure all workers are born before we quit, + // otherwise the worker will wait infinitely for work. + break; + } + } + + DEBUG_MSG(1.0, "DISPATCH: Waiting for workers to die..."); + + + /* Wait for all to die */ + while (n_slaves_alive != 0) { + // we received a message from the last loop + + if (message == GIVE_ME_WORK) { + char buf[1]; + + buf[0] = 0; // tell them to die + MPI_Send(buf, sizeof(buf), MPI_CHAR, rank, + runner_->tag_, MPI_COMM_WORLD); + + DEBUG_MSG(1.0, "DISPATCH: I told a worker to die."); + n_slaves_busy--; + } else if (message == DEATH) { + DEBUG_MSG(1.0, "DISPATCH: A worker died."); + n_slaves_alive--; + } else { + DEBUG_ASSERT_MSG(0, "DISPATCHED: Message was %d??", + message); + } + + if (n_slaves_alive != 0) { + runner_->RecvInt_(&rank, &message); + } + } + + DEBUG_MSG(1.0, "DISPATCH: All workers have died."); + + delete this; + } + }; + + private: + MPIGrainRunner *runner_; + GrainQueue *queue_; + + public: + Dispatcher() {} + ~Dispatcher() {} + + void Init(MPIGrainRunner *runner_in) { + queue_ = NULL; + runner_ = runner_in; + } + + void set_queue(GrainQueue *queue_in) { + queue_ = queue_in; + } + + private: + void MasterLoop_() { + if (runner_->n_slaves_ > 0) { + MPIMasterTask *task = new MPIMasterTask(runner_); + task->Run(); + } + } + }; + + private: + class ConsumerTask : public Task { + FORBID_COPY(ConsumerTask); + + private: + struct MPIGrainRunner *runner_; + + public: + ConsumerTask(MPIGrainRunner *runner_in) { + runner_ = runner_in; + } + + void Run() { + int my_grains = 0; + for (;;) { + Grain *grain = runner_->NextGrain_(); + if (!grain) { + break; + } + grain->Run(runner_->context_); + delete grain; + my_grains++; + } + DEBUG_MSG(1.0, "A thread on rank %d is done, completed %d grains.", + runner_->my_rank_, my_grains); + delete this; + } + }; + + private: + enum { BIRTH = 27, GIVE_ME_WORK, DEATH }; + + private: + int tag_; + Context context_; + int master_rank_; + int my_rank_; + int n_nodes_; + int n_slaves_; + Dispatcher *dispatcher_; + + WaitCondition need_work_cond_; + volatile int need_work_; + WaitCondition have_work_cond_; + volatile int have_work_; + Mutex mutex_; + ArrayList slave_grains_; + + public: + MPIGrainRunner() {} + ~MPIGrainRunner() {} + + /** + * Initialize this. + * + * @param name name to associate information with + * @param tag_in the tag for this to use; this will occupy both the + * specified tag AND the one after it + * @param context_in the context to run grains with + * @param master_rank_in the rank of the master (defaults to first) + */ + void Init(const char *name, int tag_in, Context context_in, int master_rank_in = 0) { + tag_ = tag_in; + context_ = context_in; + master_rank_ = master_rank_in; + MPI_Comm_rank(MPI_COMM_WORLD, &my_rank_); + MPI_Comm_size(MPI_COMM_WORLD, &n_nodes_); + n_slaves_ = n_nodes_ - 1; + DEBUG_MSG(1.0, "my_rank = %d", my_rank_); + if (master_rank_ == my_rank_) { + dispatcher_ = new Dispatcher(); + dispatcher_->Init(this); + } else { + dispatcher_ = NULL; + } + have_work_ = 0; + need_work_ = 0; + slave_grains_.Init(); + xrun_subparam_set(name, "n_nodes", "%d", n_nodes_); + } + + Dispatcher *dispatcher() const { + return dispatcher_; + } + + Thread *SpawnThread() { + ConsumerTask *task = new ConsumerTask(this); + Thread *thread = new Thread(); + thread->Init(task); + thread->Start(); + return thread; + } + + /** + * Creates the specified number of threads, and uses those to execute + * all grains of work. + * + * Use this rather than the other SpawnThread methods. + */ + void RunThreads(int num_threads) { + DEBUG_MSG(1.0, "Rank %d is ready to roll.", my_rank_); + ArrayList threads; + int num_worker_threads = num_threads; + + if (num_worker_threads != 2) abort(); + threads.Init(num_worker_threads); + + for (int i = 0; i < num_worker_threads; i++) { + threads[i] = SpawnThread(); + } + if (dispatcher_) { + dispatcher_->MasterLoop_(); + DEBUG_MSG(1.0, "Master loop done, cleaning up."); + } else { + SlaveLoop_(); + DEBUG_MSG(1.0, "Slave %d loop done, cleaning up.", my_rank_); + } + for (int i = num_worker_threads; i--;) { + threads[i]->WaitStop(); + delete threads[i]; + } + if (!dispatcher_) { + // Let the master know I'm done. + SendInt_(master_rank_, DEATH); + } + DEBUG_MSG(1.0, "Rank %d killed all threads.", my_rank_); + } + + private: + void SendInt_(int dest_rank, int num) { + MPI_Send(&num, 1, + MPI_INT, dest_rank, + tag_ + 1, + MPI_COMM_WORLD); + } + + void RecvInt_(int *send_rank, int *num_ptr) { + MPI_Status status; + MPI_Recv(num_ptr, 1, + MPI_INT, + MPI_ANY_SOURCE, + tag_ + 1, + MPI_COMM_WORLD, + &status); + *send_rank = status.MPI_SOURCE; + } + + Grain *NextGrain_() { + Grain *grain; + + if (dispatcher_) { + DEBUG_ASSERT((my_rank_ == master_rank_)); + mutex_.Lock(); + grain = dispatcher_->queue_->Pop(); + mutex_.Unlock(); + } else { + DEBUG_ASSERT((my_rank_ != master_rank_)); + + mutex_.Lock(); + need_work_++; + mutex_.Unlock(); + + need_work_cond_.Signal(); + + mutex_.Lock(); + while (!have_work_) { + have_work_cond_.Wait(&mutex_); + } + if (slave_grains_.size() > 0) { + grain = *slave_grains_.PopBackPtr(); + have_work_--; + } else { + grain = NULL; + } + mutex_.Unlock(); + DEBUG_MSG(2.0, "Slave gave me stuff!"); + } + + return grain; + } + + void SlaveLoop_() { + bool done = false; + + DEBUG_MSG(1.0, "%d, WORKER: Announcing birth...", my_rank_); + SendInt_(master_rank_, BIRTH); + + while (!done) { + char buf[sizeof(Grain) + 1] = "q"; + MPI_Status status; + + DEBUG_MSG(1.0, "%d, WORKER: Waiting for work...", my_rank_); + MPI_Recv(buf, sizeof(buf), + MPI_CHAR, MPI_ANY_SOURCE, + tag_, + MPI_COMM_WORLD, &status); + + Grain *grain = NULL; + + if (buf[0] == 1) { + DEBUG_MSG(1.0, "%d, WORKER: Got some work. There are %d waiting.", + my_rank_, need_work_); + + grain = new Grain(); + mem::CopyBytes(grain, buf+1, sizeof(Grain)); + + mutex_.Lock(); + while (need_work_ == 0) { + need_work_cond_.Wait(&mutex_); + } + need_work_--; + have_work_++; + *slave_grains_.AddBack() = grain; + mutex_.Unlock(); + + have_work_cond_.Signal(); + + SendInt_(master_rank_, GIVE_ME_WORK); + } else { + DEBUG_ASSERT_MSG(buf[0] == 0, "buf[0] = %d", buf[0]); + grain = NULL; + done = true; + } + } + + DEBUG_MSG(1.0, "%d, WORKER: Duly dying !!!!!!!!!!!!!!!!!!!", my_rank_); + + mutex_.Lock(); + have_work_ = -1; + mutex_.Unlock(); + + have_work_cond_.Broadcast(); + } +}; + +#endif diff --git a/fastlib/par/task.h b/fastlib/par/task.h new file mode 100644 index 0000000000..5396b6ef91 --- /dev/null +++ b/fastlib/par/task.h @@ -0,0 +1,30 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @file task.h + * + * Declaration for generic task concept. + */ + +#ifndef PAR_TASK_H +#define PAR_TASK_H + +#include "base/cc.h" +#include "base/common.h" + +/** + * Single start-to-finish task to be executed. + * + * This is a polymorphic class. + */ +class Task { + FORBID_COPY(Task); + + public: + Task() {} + virtual ~Task() {} + + virtual void Run() = 0; +}; + +#endif diff --git a/fastlib/par/thread.cc b/fastlib/par/thread.cc new file mode 100644 index 0000000000..d47c299c25 --- /dev/null +++ b/fastlib/par/thread.cc @@ -0,0 +1,5 @@ +#include "thread.h" + +Mutex Mutex::global; + +// TODO: Blank file diff --git a/fastlib/par/thread.h b/fastlib/par/thread.h new file mode 100644 index 0000000000..3848d5bfcf --- /dev/null +++ b/fastlib/par/thread.h @@ -0,0 +1,203 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @file thread.h + * + * Abstractions for helping you write threaded programs. + */ +#ifndef PAR_THREAD_H +#define PAR_THREAD_H + +#include "task.h" + +#include "base/common.h" + +#include + +/** + * Thread convenience wrapper. + * + * Usage: Create a Thread, give it a Task object, tell the thread to run, + * and eventually wait for the thread to finish. + */ +class Thread { + FORBID_COPY(Thread); + + private: +#ifdef DEBUG + enum {UNINIT, READY, ATTACHED, DETACHED, DONE} status; +#endif + pthread_t thread_; + Task *task_; + + static void *ThreadMain_(void *self) { + Thread* thread = reinterpret_cast(self); + + thread->task_->Run(); + + return NULL; + } + + void Exit_() { + pthread_exit(NULL); + } + + public: + Thread() { + DEBUG_ONLY(status = UNINIT); + } + ~Thread() { + DEBUG_ASSERT(status == DETACHED || status == READY || status == DONE); + DEBUG_ONLY(status = UNINIT); + } + + /** + * Initializes, given a task to run. + */ + void Init(Task* task_in) { + DEBUG_ASSERT(status == UNINIT); + task_ = task_in; + DEBUG_ONLY(status = READY); + } + + /** + * Starts the thread running. + */ + void Start() { + DEBUG_ASSERT(status == READY); + pthread_create(&thread_, NULL, + ThreadMain_, reinterpret_cast(this)); + DEBUG_ONLY(status = ATTACHED); + } + + /** + * Detaches a thread -- the thread will cease to exist once the task + * completes. You may not call WaitStop on this thread afterwards. + */ + void Detach() { + DEBUG_ASSERT(status == ATTACHED); + pthread_detach(thread_); + DEBUG_ONLY(status = DETACHED); + } + + /** + * Wait for a thread to stop. + * + * Failure to do this may cause your program to hang when it is done. + */ + void WaitStop() { + DEBUG_ASSERT(status == ATTACHED); + pthread_join(thread_, NULL); + DEBUG_ONLY(status = DONE); + } + + /** + * Gets the contained task. + */ + Task* task() const { + return task_; + } +}; + +/** + * Mutual exclusion lock to prevent threads from clobbering results. + */ +class Mutex { + FORBID_COPY(Mutex); + friend class WaitCondition; + + private: + pthread_mutex_t mutex_; + + public: + static Mutex global; + + public: + Mutex() { + pthread_mutex_init(&mutex_, NULL); + } + + ~Mutex() { + pthread_mutex_destroy(&mutex_); + } + + /** Obtains the lock. */ + void Lock() { + pthread_mutex_lock(&mutex_); + } + + /** Tries to lock, returns false if doing so would require waiting. */ + bool TryLock() { + return likely(!pthread_mutex_trylock(&mutex_)); + } + + /** Releases the lock. */ + void Unlock() { + pthread_mutex_unlock(&mutex_); + } +}; + +/** + * Wait condition for alerting other threads of an action. + */ +class WaitCondition { + FORBID_COPY(WaitCondition); + + private: + pthread_cond_t cond_; + + public: + WaitCondition() { + pthread_cond_init(&cond_, NULL); + } + + ~WaitCondition() { + pthread_cond_destroy(&cond_); + } + + void Signal() { + pthread_cond_signal(&cond_); + } + + void Broadcast() { + pthread_cond_broadcast(&cond_); + } + + void Wait(Mutex* mutex_to_unlock) { + pthread_cond_wait(&cond_, &mutex_to_unlock->mutex_); + } + + void WaitMillis(Mutex& mutex_to_unlock, unsigned millis) { + struct timespec ts; + + ts.tv_sec = millis / 1000; + ts.tv_nsec = (millis % 1000) * 1000000; + + pthread_cond_timedwait(&cond_, &mutex_to_unlock.mutex_, &ts); + } + + void WaitSec(Mutex& mutex_to_unlock, unsigned sec) { + struct timespec ts; + + ts.tv_sec = sec; + ts.tv_nsec = 0; + + pthread_cond_timedwait(&cond_, &mutex_to_unlock.mutex_, &ts); + } +}; + +/** + * Mix-in to make a version of an existing object that can be locked. + * + * Your object must have default constructors and use Init methods. + * The resulting object will have Lock, Unlock, and TryLock methods. + */ +template +class Lockable : public TContained, public Mutex { + FORBID_COPY(Lockable); + + Lockable() {} + ~Lockable() {} +}; + +#endif diff --git a/fastlib/tree/bounds.h b/fastlib/tree/bounds.h new file mode 100644 index 0000000000..5ddeb76b84 --- /dev/null +++ b/fastlib/tree/bounds.h @@ -0,0 +1,353 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @param bounds.h + * + * Bounds that are useful for binary space partitioning trees. + * + * TODO: Come up with a better design so you can do plug-and-play distance + * metrics. + */ + +#ifndef TREE_BOUNDS_H +#define TREE_BOUNDS_H + +#include "la/matrix.h" +#include "la/la.h" + +/** + * Simple real-valued range. + */ +struct DBound { + public: + double lo; + double hi; + + public: + DBound() {} + + void Init() { + lo = DBL_MAX; + hi = -DBL_MAX; + } + + double width() const { + return hi - lo; + } + + double mid() const { + return (hi + lo) / 2; + } +}; + +/** + * Hyper-rectangle bound. + */ +class DHrectBound { + private: + DBound *bounds_; + //double diagonal_sq_; + index_t dim_; + + public: + DHrectBound() { + DEBUG_POISON_PTR(bounds_); + DEBUG_ONLY(dim_ = BIG_BAD_NUMBER); + } + + ~DHrectBound() { + mem::Free(bounds_); + } + + template + void Deserialize(Deserializer *s) { + DEBUG_ASSERT_MSG(dim_ == BIG_BAD_NUMBER, "Already initialized"); + + s->Get(&dim_); + bounds_ = mem::Alloc(dim_); + s->Get(bounds_, dim_); + + //ComputeDiagonal_(); + } + + template + void Serialize(Serializer *s) const { + s->Put(dim_); + s->Put(bounds_, dim_); + } + + void Init(index_t dimension) { + DEBUG_ASSERT_MSG(dim_ == BIG_BAD_NUMBER, "Already initialized"); + + bounds_ = mem::Alloc(dimension); + + for (index_t i = 0; i < dimension; i++) { + bounds_[i].Init(); + } + + dim_ = dimension; + + //ComputeDiagonal_(); + } + + bool Belongs(const Vector& point) const { + for (index_t i = 0; i < point.length(); i++) { + const DBound *bound = &bounds_[i]; + if (point[i] > bound->hi || point[i] < bound->lo) { + return false; + } + } + + return true; + } + + double MinDistanceSqToInstance(const Vector& point) const { + DEBUG_ASSERT(point.length() == dim_); + return MinDistanceSqToInstance(point.ptr()); + } + + double MinDistanceSqToInstance(const double *mpoint) const { + double sumsq = 0; + //index_t mdim = dim_; + const DBound *mbound = bounds_; + + index_t d = dim_; + + do { + double v = *mpoint; + double v1 = mbound->lo - v; + double v2 = v - mbound->hi; + + v = (v1 + fabs(v1)) + (v2 + fabs(v2)); + + mbound++; + mpoint++; + + sumsq += v * v; + } while (--d); + + return sumsq / 4; + } + + double MaxDistanceSqToInstance(const Vector& point) const { + double sumsq = 0; + + DEBUG_ASSERT(point.length() == dim_); + + for (index_t d = 0; d < dim_; d++) { + double v = max(point[d] - bounds_[d].lo, + bounds_[d].hi - point[d]); + + sumsq += v * v; + } + + return sumsq; + } + + double MinDistanceSqToBound(const DHrectBound& other) const { + double sumsq = 0; + const DBound *a = this->bounds_; + const DBound *b = other.bounds_; + index_t mdim = dim_; + + DEBUG_ASSERT(dim_ == other.dim_); + + // We invoke the following: + // x + fabs(x) = max(x * 2, 0) + // (x * 2)^2 / 4 = x^2 + + for (index_t d = 0; d < mdim; d++) { +#if 0 + double v = b[d].lo - a[d].hi; + + if (v < 0) { + v = a[d].lo - b[d].hi; + } + + if (likely(v > 0)) { + sumsq += v * v; + } +#else + double v1 = b[d].lo - a[d].hi; + double v2 = a[d].lo - b[d].hi; + + double v = (v1 + fabs(v1)) + (v2 + fabs(v2)); + + sumsq += v * v; +#endif + } + + return sumsq / 4; + } + + double MinDistanceSqToBoundFarEnd(const DHrectBound& other) const { + double sumsq = 0; + const DBound *a = this->bounds_; + const DBound *b = other.bounds_; + index_t mdim = dim_; + + DEBUG_ASSERT(dim_ == other.dim_); + + for (index_t d = 0; d < mdim; d++) { + double v1 = b[d].hi - a[d].hi; + double v2 = a[d].lo - b[d].lo; + + double v = max(v1, v2); + v = (v + fabs(v)); /* truncate negative */ + + sumsq += v * v; + } + + return sumsq / 4; + } + + double MaxDistanceSqToBound(const DHrectBound& other) const { + double sumsq = 0; + const DBound *a = this->bounds_; + const DBound *b = other.bounds_; + + DEBUG_ASSERT(dim_ == other.dim_); + + for (index_t d = 0; d < dim_; d++) { + double v = max(b[d].hi - a[d].lo, a[d].hi - b[d].lo); + + sumsq += v * v; + } + + return sumsq; + } + + double MidDistanceSqToBound(const DHrectBound& other) const { + double sumsq = 0; + const DBound *a = this->bounds_; + const DBound *b = other.bounds_; + + DEBUG_ASSERT(dim_ == other.dim_); + + for (index_t d = 0; d < dim_; d++) { + double v = (a[d].hi + a[d].lo - b[d].hi - b[d].lo) * 0.5; + + sumsq += v * v; + } + + return sumsq; + } + + void Update(const Vector& vector) { + DEBUG_ASSERT(vector.length() == dim_); + + for (index_t i = 0; i < dim_; i++) { + DBound* bound = &bounds_[i]; + double d = vector[i]; + + if (unlikely(d > bound->hi)) { + bound->hi = d; + } + if (unlikely(d < bound->lo)) { + bound->lo = d; + } + } + } + + const DBound& get(index_t i) const { + return bounds_[i]; + } + + //double diagonal_sq() const { + // return diagonal_sq_; + //} + + FORBID_COPY(DHrectBound); + + private: + //void ComputeDiagonal_() { + // diagonal_sq_ = 0; + // for (index_t d = 0; d < dim_; d++) { + // double v = bounds_[d].lo - bounds_[d].hi; + // diagonal_sq_ += v*v; + // } + //} +}; + +/** + * Euclidean metric for use with ball bounds. + */ +class DEuclideanMetric { + public: + static double CalculateMetric(const Vector& a, const Vector& b) { + return sqrt(la::DistanceSqEuclidean(a.ptr(), b.ptr(), a.length())); + } +}; + +/** + * Bound of a ball tree. + */ +template +class BallBound { + FORBID_COPY(BallBound); + + public: + typedef TMetric Metric; + typedef TInstance Instance; + + private: + Instance center_; + double radius_; + + public: + BallBound() {} + + const Instance& center() const { + return center; + } + + Instance& center() { + return center; + } + + double radius() const { + return radius; + } + + void set_radius(double d) { + radius = d; + } + + double DistanceToCenter(const Instance& point) { + return Metric::CalculateMetric(point, center_); + } + + bool Belongs(const Instance& point) { + return DistanceToCenter(point) <= radius_; + } + + double MinDistanceToInstance(const Instance& point) { + return max(0.0, DistanceToCenter(point) - radius_); + } + + double MaxDistanceToInstance(const Instance& point) { + return DistanceToCenter(point) + radius_; + } + + double MinDistanceToBound(const BallBound& ball) { + return max(0, + DistanceToCenter(ball.center_) - (radius_ + ball.radius_)); + } + + double MaxDistanceToBound(const BallBound& ball) { + return DistanceToCenter(ball.center_) + (radius_ + ball.radius_); + } + + double MidDistanceToBound(const BallBound& other) { + return DistanceToCenter(other.center_); + } + + double MidDistanceToInstance(const Instance& point) { + return DistanceToCenter(point); + } +}; + +typedef BallBound DEuclideanBallBound; + + +#endif diff --git a/fastlib/tree/build.py b/fastlib/tree/build.py new file mode 100644 index 0000000000..ca8eed872a --- /dev/null +++ b/fastlib/tree/build.py @@ -0,0 +1,10 @@ + +librule( + sources = [], + headers = [ + "kdtree.h", "bounds.h", "spacetree.h", "statistic.h" + ], + deplibs = ["base:base", "la:la", "col:col", + "file:file"] + ) + diff --git a/fastlib/tree/kdtree.h b/fastlib/tree/kdtree.h new file mode 100644 index 0000000000..9a843bfbbc --- /dev/null +++ b/fastlib/tree/kdtree.h @@ -0,0 +1,253 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @file kdtree.h + * + * Tools for kd-trees. + * + * Eventually we hope to support KD trees with non-L2 (Euclidean) + * metrics, like Manhattan distance. + */ + +#ifndef TREE_KDTREE_H +#define TREE_KDTREE_H + +#include "spacetree.h" +#include "bounds.h" + +#include "base/common.h" +#include "col/arraylist.h" +#include "file/serialize.h" + +/* Implementation */ + +namespace tree_kdtree_private { + + template + void FindBoundFromMatrix(const Matrix& matrix, + index_t first, index_t count, TBound *bounds) { + index_t end = first + count; + for (index_t i = first; i < end; i++) { + Vector col; + + matrix.MakeColumnVector(i, &col); + bounds->Update(col); + } + } + + template + index_t MatrixPartition( + Matrix& matrix, index_t dim, double splitvalue, + index_t first, index_t count, + TBound* left_bound, TBound* right_bound, + index_t *old_from_new) { + index_t left = first; + index_t right = first + count - 1; + + /* At any point: + * + * everything < left is correct + * everything > right is correct + */ + for (;;) { + while (matrix.get(dim, left) < splitvalue && likely(left <= right)) { + Vector left_vector; + matrix.MakeColumnVector(left, &left_vector); + left_bound->Update(left_vector); + left++; + } + + while (matrix.get(dim, right) >= splitvalue && likely(left <= right)) { + Vector right_vector; + matrix.MakeColumnVector(right, &right_vector); + right_bound->Update(right_vector); + right--; + } + + if (unlikely(left > right)) { + /* left == right + 1 */ + break; + } + + Vector left_vector; + Vector right_vector; + + matrix.MakeColumnVector(left, &left_vector); + matrix.MakeColumnVector(right, &right_vector); + + left_vector.SwapValues(&right_vector); + + left_bound->Update(left_vector); + right_bound->Update(right_vector); + + if (old_from_new) { + index_t t = old_from_new[left]; + old_from_new[left] = old_from_new[right]; + old_from_new[right] = t; + } + + DEBUG_ASSERT(left <= right); + right--; + + // this conditional is always true, I belueve + //if (likely(left <= right)) { + // right--; + //} + } + + DEBUG_ASSERT(left == right + 1); + + return left; + } + + template + void SplitKdTreeMidpoint(Matrix& matrix, + TKdTree *node, index_t leaf_size, index_t *old_from_new) { + TKdTree *left = NULL; + TKdTree *right = NULL; + + //FindBoundFromMatrix(matrix, node->begin(), node->count(), + // &node->bound()); + + if (node->count() > leaf_size) { + index_t split_dim = BIG_BAD_NUMBER; + double max_width = -1; + + for (index_t d = 0; d < matrix.n_rows(); d++) { + double w = node->bound().get(d).width(); + + if (unlikely(w > max_width)) { + max_width = w; + split_dim = d; + } + } + + double split_val = node->bound().get(split_dim).mid(); + + if (max_width == 0) { + // Okay, we can't do any splitting, because all these points are the + // same. We have to give up. + } else { + left = new TKdTree(); + left->bound().Init(matrix.n_rows()); + + right = new TKdTree(); + right->bound().Init(matrix.n_rows()); + + index_t split_col = MatrixPartition(matrix, split_dim, split_val, + node->begin(), node->count(), + &left->bound(), &right->bound(), + old_from_new); + + DEBUG_MSG(3.0,"split (%d,[%d],%d) dim %d on %f (between %f, %f)", + node->begin(), split_col, + node->begin() + node->count(), split_dim, split_val, + node->bound().get(split_dim).lo, + node->bound().get(split_dim).hi); + + left->Init(node->begin(), split_col - node->begin()); + right->Init(split_col, node->begin() + node->count() - split_col); + + // This should never happen if max_width > 0 + DEBUG_ASSERT(left->count() != 0 && right->count() != 0); + + SplitKdTreeMidpoint(matrix, left, leaf_size, old_from_new); + SplitKdTreeMidpoint(matrix, right, leaf_size, old_from_new); + } + } + + node->set_children(matrix, left, right); + } +}; + +namespace tree { + /** + * Creates a KD tree from data, splitting on the midpoint. + * + * This requires you to pass in two unitialized ArrayLists which will contain + * index mappings so you can account for the re-ordering of the matrix. + * (By unitialized I mean don't call Init on it) + * + * @param matrix data where each column is a point, WHICH WILL BE RE-ORDERED + * @param old_from_new pointer to an unitialized arraylist; it will map + * original indexes to new indices + * @param old_from_new pointer to an unitialized arraylist; it will map + * new indices to original + */ + template + TKdTree *MakeKdTreeMidpoint(Matrix& matrix, index_t leaf_size, + ArrayList *old_from_new = NULL, + ArrayList *new_from_old = NULL) { + TKdTree *node = new TKdTree(); + index_t *old_from_new_ptr; + + if (old_from_new) { + old_from_new->Init(matrix.n_cols()); + + for (index_t i = 0; i < matrix.n_cols(); i++) { + (*old_from_new)[i] = i; + } + + old_from_new_ptr = old_from_new->begin(); + } else { + old_from_new_ptr = NULL; + } + + node->Init(0, matrix.n_cols()); + node->bound().Init(matrix.n_rows()); + tree_kdtree_private::FindBoundFromMatrix(matrix, + 0, matrix.n_cols(), &node->bound()); + + tree_kdtree_private::SplitKdTreeMidpoint(matrix, node, leaf_size, + old_from_new_ptr); + + if (new_from_old) { + new_from_old->Init(matrix.n_cols()); + for (index_t i = 0; i < matrix.n_cols(); i++) { + (*new_from_old)[(*old_from_new)[i]] = i; + } + } + + return node; + } + + // TODO: Perhaps move this into a "util.h" file + + template + void SerializeKdTree(const TKdTree *tree, + const Matrix& matrix, + const ArrayList& old_from_new, + Serializer *s) { + s->PutMagic(file::CreateMagic("kdtree")); + tree->SerializeAll(matrix, s); + old_from_new.Serialize(s); + } + + template + void DeserializeKdTree(TKdTree *uninit_tree, + Matrix* uninit_matrix, + ArrayList* uninit_old_from_new, + Deserializer *s) { + s->AssertMagic(file::CreateMagic("kdtree")); + uninit_tree->DeserializeAll(uninit_matrix, s); + if (uninit_old_from_new) { + uninit_old_from_new->Deserialize(s); + } + } + + template + void ReadKdTreeFromFile(const char *fname, + TKdTree *uninit_tree, + Matrix* uninit_matrix, + ArrayList* uninit_old_from_new) { + NativeFileDeserializer ds; + + ASSERT_PASS(ds.Init(fname)); + DeserializeKdTree(uninit_tree, uninit_matrix, uninit_old_from_new, + &ds); + } +}; + +typedef BinarySpaceTree BasicKdTree; + +#endif diff --git a/fastlib/tree/spacetree.h b/fastlib/tree/spacetree.h new file mode 100644 index 0000000000..5d8c01a686 --- /dev/null +++ b/fastlib/tree/spacetree.h @@ -0,0 +1,279 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @file spacetree.h + * + * Generalized space partitioning tree. + */ + +#ifndef TREE_SPACETREE_H +#define TREE_SPACETREE_H + +#include "base/cc.h" +#include "file/serialize.h" +#include "statistic.h" + +/** + * A binary space partitioning tree, such as KD or ball tree. + * + * This particular tree forbids you from having more children. + * + * @param TBound the bounding type of each child (TODO explain interface) + * @param TDataset the data set type + * @param TStatistic extra data in the node + */ +template > +class BinarySpaceTree { + public: + typedef TBound Bound; + typedef TDataset Dataset; + typedef TStatistic Statistic; + + private: + Bound bound_; + BinarySpaceTree *left_; + BinarySpaceTree *right_; + index_t begin_; + index_t count_; + Statistic stat_; + + public: + BinarySpaceTree() { + DEBUG_ONLY(begin_ = BIG_BAD_NUMBER); + DEBUG_ONLY(count_ = BIG_BAD_NUMBER); + DEBUG_POISON_PTR(left_); + DEBUG_POISON_PTR(right_); + } + + + ~BinarySpaceTree() { + if (!is_leaf()) { + delete left_; + delete right_; + } + DEBUG_ONLY(begin_ = BIG_BAD_NUMBER); + DEBUG_ONLY(count_ = BIG_BAD_NUMBER); + DEBUG_POISON_PTR(left_); + DEBUG_POISON_PTR(right_); + } + + void Init(index_t begin_in, index_t count_in) { + DEBUG_ASSERT(begin_ == BIG_BAD_NUMBER); + DEBUG_POISON_PTR(left_); + DEBUG_POISON_PTR(right_); + begin_ = begin_in; + count_ = count_in; + } + + /** + * Find a node in this tree by its begin and count. + * + * Every node is uniquely identified by these two numbers. + * This is useful for communicating position over the network, + * when pointers would be invalid. + * + * @param begin_q the begin() of the node to find + * @param count_q the count() of the node to find + * @return the found node, or NULL + */ + const BinarySpaceTree* FindByBeginCount( + index_t begin_q, index_t count_q) const { + DEBUG_ASSERT(begin_q >= begin_); + DEBUG_ASSERT(count_q <= count_); + if (begin_ == begin_q && count_ == count_q) { + return this; + } else if (unlikely(is_leaf())) { + return NULL; + } else if (begin_q < right_->begin_) { + return left_->FindByBeginCount(begin_q, count_q); + } else { + return right_->FindByBeginCount(begin_q, count_q); + } + } + + /** + * Find a node in this tree by its begin and count (const). + * + * Every node is uniquely identified by these two numbers. + * This is useful for communicating position over the network, + * when pointers would be invalid. + * + * @param begin_q the begin() of the node to find + * @param count_q the count() of the node to find + * @return the found node, or NULL + */ + BinarySpaceTree* FindByBeginCount( + index_t begin_q, index_t count_q) { + DEBUG_ASSERT(begin_q >= begin_); + DEBUG_ASSERT(count_q <= count_); + if (begin_ == begin_q && count_ == count_q) { + return this; + } else if (unlikely(is_leaf())) { + return NULL; + } else if (begin_q < right_->begin_) { + return left_->FindByBeginCount(begin_q, count_q); + } else { + return right_->FindByBeginCount(begin_q, count_q); + } + } + + /** + * Serializes the tree _structure_ only. + * Statistics are not stored (this allows you to re-load the tree + * for problems that require different statistics). + */ + template + void Serialize(Serializer *s) const { + bound_.Serialize(s); + s->Put(begin_); + s->Put(count_); + + bool children = !is_leaf(); + s->Put(children); + + if (children) { + left_->Serialize(s); + right_->Serialize(s); + } + } + + /** + * Deserializes the tree from its structure, and re-calculates bottom-up + * statistics. + */ + template + void Deserialize(const Dataset& data, Deserializer *s) { + DEBUG_ASSERT(begin_ == BIG_BAD_NUMBER); + + bound_.Deserialize(s); + s->Get(&begin_); + s->Get(&count_); + + bool children; + + s->Get(&children); + + BinarySpaceTree *l = NULL; + BinarySpaceTree *r = NULL; + + if (children) { + l = new BinarySpaceTree(); + l->Deserialize(data, s); + r = new BinarySpaceTree(); + r->Deserialize(data, s); + } + + set_children(data, l, r); + } + + template + void SerializeAll(const Dataset& data, Serializer *s) const { + // can't use BinarySpaceTree as a magic number, because we want to be + // able to deserialize this class with another statistic + s->PutMagic(file::CreateMagic("spacetree") + + MAGIC_NUMBER(TDataset) + MAGIC_NUMBER(TBound)); + data.Serialize(s); + Serialize(s); + } + + template + void DeserializeAll(Dataset* data, Deserializer *s) { + // can't use BinarySpaceTree as a magic number, because we want to be + // able to deserialize this class with another statistic + s->AssertMagic(file::CreateMagic("spacetree") + + MAGIC_NUMBER(TDataset) + MAGIC_NUMBER(TBound)); + data->Deserialize(s); + Deserialize(*data, s); + } + + // TODO: Not const correct + + /** + * Used only when constructing the tree. + */ + void set_children(const Dataset& data, + BinarySpaceTree *left_in, BinarySpaceTree *right_in) { + left_ = left_in; + right_ = right_in; + if (!is_leaf()) { + stat_.Init(data, begin_, count_, left_->stat_, right_->stat_); + DEBUG_ASSERT(count_ == left_->count_ + right_->count_); + DEBUG_ASSERT(left_->begin_ == begin_); + DEBUG_ASSERT(right_->begin_ == begin_ + left_->count_); + } else { + stat_.Init(data, begin_, count_); + } + } + + const Bound& bound() const { + return bound_; + } + + Bound& bound() { + return bound_; + } + + const Statistic& stat() const { + return stat_; + } + + Statistic& stat() { + return stat_; + } + + bool is_leaf() const { + return !left_; + } + + /** + * Gets the left branch of the tree. + */ + BinarySpaceTree *left() const { + // TODO: Const correctness + return left_; + } + + /** + * Gets the right branch. + */ + BinarySpaceTree *right() const { + // TODO: Const correctness + return right_; + } + + /** + * Gets the index of the begin point of this subset. + */ + index_t begin() const { + return begin_; + } + + /** + * Gets the index one beyond the last index in the series. + */ + index_t end() const { + return begin_ + count_; + } + + /** + * Gets the number of points in this subset. + */ + index_t count() const { + return count_; + } + + void Print() const { + printf("node: %d to %d: %d points total\n", + begin_, begin_ + count_ - 1, count_); + if (!is_leaf()) { + left_->Print(); + right_->Print(); + } + } + + FORBID_COPY(BinarySpaceTree); +}; + +#endif diff --git a/fastlib/tree/statistic.h b/fastlib/tree/statistic.h new file mode 100644 index 0000000000..c6d1b64b36 --- /dev/null +++ b/fastlib/tree/statistic.h @@ -0,0 +1,40 @@ +// Copyright 2007 Georgia Institute of Technology. All rights reserved. +// ABSOLUTELY NOT FOR DISTRIBUTION +/** + * @file statistic.h + * + * Home for the concept of tree statistics. + * + * You should define your own statistic that looks like EmptyStatistic. + */ + +#ifndef TREE_STATISTIC_H +#define TREE_STATISTIC_H + +/** + * Empty statistic if you are not interested in storing statistics in your + * tree. Use this as a template for your own. + */ +template +class EmptyStatistic { + public: + EmptyStatistic() {} + ~EmptyStatistic() {} + + /** + * Initializes by taking statistics on raw data. + */ + void Init(const TDataset& dataset, index_t start, index_t count) { + } + + /** + * Initializes by combining statistics of two partitions. + * + * This lets you build fast bottom-up statistics when building trees. + */ + void Init(const TDataset& dataset, index_t start, index_t count, + const EmptyStatistic& left_stat, const EmptyStatistic& right_stat) { + } +}; + +#endif