Initial Branch Commit

list.c

+#include "list.h"
+
+
+const int GROWTH_F = 2;
+
+typedef struct list
+{
+    void **data;
+    size_t size;
+    size_t cap;
+    free_func_t free_func;
+    pthread_mutex_t lock;
+    pthread_cond_t cond;
+} list_t;
+
+list_t *list_init(size_t initial_cap, free_func_t freer)
+{
+    list_t *list = malloc(sizeof(list_t));
+    pthread_mutex_init(&list->lock, NULL);
+    pthread_cond_init(&list->cond, NULL);
+    if (initial_cap < 1)
+    {
+        list->cap = 1;
+    }
+    else
+    {
+        list->cap = initial_cap;
+    }
+
+    list->size = 0;
+    list->data = malloc(sizeof(void *) * list->cap);
+    list->free_func = freer;
+    return list;
+}
+
+void list_resize(list_t *list)
+{
+  list->data = realloc(list->data, sizeof(void *) * GROWTH_F * list->cap);
+  assert(list->data != NULL);
+  list->cap = list->cap * GROWTH_F;
+}
+
+void list_free(list_t *list)
+{
+    pthread_mutex_lock(&list->lock);
+    if (list->free_func != NULL)
+    {
+        for (int i = 0; i < list->size; i++)
+        {
+            if (list->data[i] != NULL)
+            {
+                (*list->free_func)(list->data[i]);
+            }
+        }
+    }
+    while (list->size >= 0) {
+        list_remove_front(list);
+    }
+    pthread_mutex_unlock(&list->lock);
+    free(list->data);
+    free(list);
+}
+
+size_t list_size(list_t *list) { return list->size; }
+
+void *list_get(list_t *list, size_t index)
+{
+    pthread_mutex_lock(&list->lock);
+    assert(index >= 0);
+    assert(index < list->size);
+    pthread_cond_signal(&list->cond);
+    pthread_mutex_unlock(&list->lock);
+    return list->data[index];
+}
+
+void list_set(list_t *list, void *value, size_t index)
+{
+    pthread_mutex_lock(&list->lock);
+    assert(list->size >= 1);
+    assert(index < list->size);
+    list->data[index] = value;
+    pthread_cond_signal(&list->cond);
+    pthread_mutex_unlock(&list->lock);
+}
+
+void *list_remove(list_t *list, size_t index)
+{
+    pthread_mutex_lock(&list->lock);
+    assert(list->size >= 1);
+    assert(index < list->size);
+    void *ret = list->data[index];
+    for (int i = index + 1; i < list->size; i++)
+    {
+        list->data[i - 1] = list->data[i];
+    }
+    list->size--;
+    pthread_cond_signal(&list->cond);
+    pthread_mutex_unlock(&list->lock);
+    return ret;
+}
+
+void *list_remove_back(list_t *list)
+{
+  return list_remove(list, list->size - 1);
+}
+
+void list_add(list_t *list, void *value)
+{
+    pthread_mutex_lock(&list->lock);
+    assert(value != NULL);
+    if (list->size == list->cap)
+    {
+        list_resize(list);
+    }
+    list->data[list->size] = value;
+    list->size++;
+    pthread_cond_signal(&list->cond);
+    pthread_mutex_unlock(&list->lock);
+}
+
+
+void *list_remove_front(list_t *list) { return list_remove(list, 0); }
+
+free_func_t list_get_freer(list_t *list) { return list->free_func; }
+
+bool list_contains(list_t *list, void *object)
+{
+  for (size_t i = 0; i > list->size; i++)
+  {
+    if (object == list->data[i])
+    {
+      return true;
+    }
+  }
+  return false;
+}
+
+list_t *list_copy(list_t *list)
+{
+  list_t *copy = list_init(list_size(list), list_get_freer(list));
+  for (size_t i = 0; i < list_size(list); i++)
+  {
+    list_add(copy, list_get(list, i));
+  }
+  return copy;
+}
+
+void **list_get_data(list_t *list){
+    return list->data;
+}

list.h

+#ifndef __LIST_H__
+#define __LIST_H__
+
+#include <stdbool.h>
+#include <stddef.h>
+#include <assert.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <pthread.h>
+/**
+ * A growable array of pointers.
+ * Can store values of any pointer type (e.g. vector_t*, body_t*).
+ * The list automatically grows its internal array when more capacity is needed.
+ */
+
+typedef struct list list_t;
+/**
+ * A function that can be called on list elements to release their resources.
+ * Examples: free, body_free
+ */
+typedef void (*free_func_t)(void *);
+
+/**
+ * Allocates memory for a new list with space for the given number of elements.
+ * The list is initially empty.
+ * Asserts that the required memory was allocated.
+ *
+ * @param initial_size the number of elements to allocate space for
+ * @param freer if non-NULL, a function to call on elements in the list
+ *   in list_free() when they are no longer in use
+ * @return a pointer to the newly allocated list
+ */
+list_t *list_init(size_t initial_size, free_func_t freer);
+
+/**
+ * Releases the memory allocated for a list.
+ *
+ * @param list a pointer to a list returned from list_init()
+ */
+void list_free(list_t *list);
+
+/**
+ * Gets the size of a list (the number of occupied elements).
+ * Note that this is NOT the list's capacity.
+ *
+ * @param list a pointer to a list returned from list_init()
+ * @return the number of elements in the list
+ */
+size_t list_size(list_t *list);
+
+/**
+ * Gets the element at a given index in a list.
+ * Asserts that the index is valid, given the list's current size.
+ *
+ * @param list a pointer to a list returned from list_init()
+ * @param index an index in the list (the first element is at 0)
+ * @return the element at the given index, as a void*
+ */
+void *list_get(list_t *list, size_t index);
+
+
+/**
+ * @brief Sets the element at the given index of the list to the value.
+ * 
+ * @param list pointer to a list returned from list_init()
+ * @param value the element to replace the existing element at the index
+ * @param index the index of the element that is being chahnged to the given value
+ */
+void list_set(list_t *list, void *value, size_t index);
+
+/**
+ * Removes the element at a given index in a list and returns it,
+ * moving all subsequent elements towards the start of the list.
+ * Asserts that the index is valid, given the list's current size.
+ *
+ * @param list a pointer to a list returned from list_init()
+ * @return the element at the given index in the list
+ */
+void *list_remove(list_t *list, size_t index);
+
+/**
+ * Removes the element at the end of a list and returns it,
+ * moving all subsequent elements towards the start of the list.
+ *
+ * @param list a pointer to a list returned from list_init()
+ * @return the element at the end of the list
+ */
+void *list_remove_back(list_t *list);
+
+/**
+ * Appends an element to the end of a list.
+ * If the list is filled to capacity, resizes the list to fit more elements
+ * and asserts that the resize succeeded.
+ * Also asserts that the value being added is non-NULL.
+ *
+ * @param list a pointer to a list returned from list_init()
+ * @param value the element to add to the end of the list
+ */
+void list_add(list_t *list, void *value);
+
+/**
+ * Removes the element at the front of a list and returns it,
+ * moving all subsequent elements towards the end of the list.
+ *
+ * @param list a pointer to a list returned from list_init()
+ * @return the element at the front of the list
+ */
+void *list_remove_front(list_t *list);
+
+/**
+ * Returns the free function passed into the list_t
+ * when initiallized.
+ *
+ * @param list a pointer to a list returned from list_init()
+ * @return the free function of the list as a free_func_t.
+ */
+free_func_t list_get_freer(list_t *list);
+
+bool list_contains(list_t *list, void *object);
+
+/**
+ * @brief
+ *
+ */
+list_t *list_copy(list_t *list);
+
+void **list_get_data(list_t *list);
+#endif // #ifndef __LIST_H__

npy.c

+// C code to write .npy files
+
+#include "npy.h"
+
+
+static void dtype_info(const int dtype, int *bytes, char *letter)
+{
+    switch (dtype)
+    {
+#define CASE(T, dtype, let) \
+    case dtype:             \
+        *bytes = sizeof(T); \
+        *letter = let;      \
+        break;
+        CASE(_Bool, NPY_BOOL, 'b')
+        CASE(signed char, NPY_BYTE, 'i')
+        CASE(unsigned char, NPY_UBYTE, 'u')
+        CASE(short, NPY_SHORT, 'i')
+        CASE(unsigned short, NPY_USHORT, 'u')
+        CASE(int, NPY_INT, 'i')
+        CASE(unsigned int, NPY_UINT, 'u')
+        CASE(long, NPY_LONG, 'i')
+        CASE(unsigned long, NPY_ULONG, 'u')
+        CASE(long long, NPY_LONGLONG, 'i')
+        CASE(unsigned long long, NPY_ULONGLONG, 'u')
+        CASE(float, NPY_FLOAT, 'f')
+        CASE(double, NPY_DOUBLE, 'f')
+        CASE(long double, NPY_LONGDOUBLE, 'f')
+#undef CASE
+    default:
+        fprintf(stderr, "npy.c: invalid dtype %d\n", dtype);
+        exit(1);
+    }
+}
+
+// Padded header size to enable easy skipping
+#define NPY_HEADER_SIZE 256
+
+void skip_npy_header(FILE *f)
+{
+    fseek(f, NPY_HEADER_SIZE, SEEK_SET);
+}
+
+void write_npy_header(FILE *f, const int ndim, const int *shape, const int dtype)
+{
+    // Grab dtype info
+    int bytes;
+    char letter;
+    dtype_info(dtype, &bytes, &letter);
+
+    // Move to beginning of file in case we just wrote the binary data
+    fseek(f, 0, SEEK_SET);
+
+    // Write header
+    const uint16_t header_size = NPY_HEADER_SIZE - 10;
+    fwrite("\x93NUMPY\x01\x00", 1, 8, f);
+    fwrite(&header_size, 1, 2, f);
+    fprintf(f, "{'descr': '<%c%d', 'fortran_order': False, 'shape': (", letter, bytes);
+    for (int i = 0; i < ndim; i++)
+        fprintf(f, "%d,", shape[i]);
+    fputs("), }", f);
+
+    // Pad with spaces out to NPY_HEADER_SIZE-1, and finish with a newline
+    fprintf(f, "%*s\n", NPY_HEADER_SIZE - 1 - (int)ftell(f), "");
+}
\ No newline at end of file

npy.h

+// C code to write .npy files
+#ifndef NPY_H 
+#define NPY_H
+
+// Example usage:
+//   FILE* f = fopen(path, "wb");
+//   skip_npy_header(f);
+//   ...write binary data...
+//   write_npy_header(f, ndim, shape, dtype);
+//   fclose(f);
+
+#include <stdio.h>
+#include <stdint.h>
+#include <stdlib.h>
+// Lifted from numpy/ndarraytypes.h
+enum NPY_TYPES
+{
+    NPY_BOOL = 0,
+    NPY_BYTE,
+    NPY_UBYTE,
+    NPY_SHORT,
+    NPY_USHORT,
+    NPY_INT,
+    NPY_UINT,
+    NPY_LONG,
+    NPY_ULONG,
+    NPY_LONGLONG,
+    NPY_ULONGLONG,
+    NPY_FLOAT,
+    NPY_DOUBLE,
+    NPY_LONGDOUBLE,
+    NPY_CFLOAT,
+    NPY_CDOUBLE,
+    NPY_CLONGDOUBLE,
+    NPY_OBJECT = 17,
+    NPY_STRING,
+    NPY_UNICODE,
+    NPY_VOID
+};
+
+// Skip space for the header at the beginning of a file.
+// Call this if you don't know the shape prior to writing the data.
+void skip_npy_header(FILE *f);
+
+// Write an .npy header at the beginning of the file
+void write_npy_header(FILE *f, const int ndim, const int *shape, const int dtype);
+
+#endif
\ No newline at end of file

queue.c

+#include "queue.h"
+
+typedef struct node {
+    void *data;
+    void *next;
+} node_t;
+
+typedef struct queue {
+    node_t *head;
+    node_t *tail;
+    pthread_mutex_t lock;
+    pthread_cond_t cond;
+    size_t size;
+    bool ready;
+} queue_t;
+
+void queue_lock(queue_t *queue) {
+    pthread_mutex_lock(&queue->lock);
+}
+
+void queue_unlock(queue_t *queue) {
+    pthread_mutex_unlock(&queue->lock);
+}
+
+void queue_sleep(queue_t *queue, bool condition) {
+    if (condition) {
+        pthread_cond_wait(&queue->cond, &queue->lock);
+    }
+}
+
+void queue_wake(queue_t *queue, bool condition) {
+    if (condition) {
+        pthread_cond_signal(&queue->cond);
+    }
+}
+
+bool queue_ready(queue_t *queue) {
+    return queue->ready;
+}
+
+node_t *node_init(void *data) {
+    node_t *node = malloc(sizeof(node_t));
+    *node = (node_t){data, NULL};
+    return node;
+}
+
+void node_free(node_t *node) {
+    if (node != NULL) {
+        free(node);
+    }
+}
+
+queue_t *queue_init(void) {
+    queue_t *queue = malloc(sizeof(queue_t));
+    queue->head = NULL;
+    queue->tail = NULL;
+    pthread_mutex_init(&queue->lock, NULL);
+    pthread_cond_init(&queue->cond, NULL);
+    queue->size = 0;
+    queue->ready = true;
+    return queue;
+}
+
+void queue_enqueue(queue_t *queue, void *value) {
+    pthread_mutex_lock(&queue->lock);
+    queue->ready = false;
+    queue->size++;
+    node_t *node = node_init(value);
+    if (queue->head == NULL && queue->tail == NULL) {
+        queue->head = node;
+        queue->tail = node;
+    }
+    else {
+        queue->tail->next = node;
+        queue->tail = node;
+    }
+    queue->ready = true;
+    pthread_cond_signal(&queue->cond);
+    pthread_mutex_unlock(&queue->lock);
+}
+
+void *queue_dequeue(queue_t *queue) {
+    pthread_mutex_lock(&queue->lock);
+    if (!queue->ready) {
+        pthread_cond_wait(&queue->cond, &queue->lock);
+    }
+    queue->ready = false;
+    while (queue->head == NULL) {
+        queue->ready = true;
+        pthread_cond_signal(&queue->cond);
+        pthread_cond_wait(&queue->cond, &queue->lock);
+        queue->ready = false;
+    }
+    node_t *head = queue->head;
+    void *out = head->data;
+    if (queue->head->next == NULL) {
+        queue->head = NULL;
+        queue->tail = NULL;
+    }
+    else {
+        queue->head = head->next;
+    }
+    queue->size--;
+    node_free(head);
+    queue->ready = true;
+    pthread_cond_signal(&queue->cond);
+    pthread_mutex_unlock(&queue->lock);
+    return out;
+}
+
+void queue_free(queue_t *queue) {
+    pthread_mutex_lock(&queue->lock);
+    queue->ready = true;
+    while (queue->head != NULL && queue->tail != NULL) {
+        queue_dequeue(queue);
+    }
+    pthread_mutex_unlock(&queue->lock);
+    free(queue);
+}

queue.h

+#ifndef QUEUE_H
+#define QUEUE_H
+
+#include <pthread.h>
+#include <stdbool.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+/** A FIFO queue */
+typedef struct queue queue_t;
+
+/**
+ * Creates a new heap-allocated FIFO queue. The queue is initially empty.
+ *
+ * @return a pointer to the new queue
+ */
+queue_t *queue_init(void);
+
+/**
+ * Enqueues a value into a queue. There is no maximum capacity,
+ * so this should succeed unless the program runs out of memory.
+ * This function should be concurrency-safe:
+ * multiple threads may call queue_enqueue() or queue_dequeue() simultaneously.
+ *
+ * @param queue the queue to append to
+ * @param value the value to add to the back of the queue
+ */
+void queue_enqueue(queue_t *queue, void *value);
+
+/**
+ * Dequeues a value from a queue.
+ * The value returned is the first enqueued value that was not yet dequeued.
+ * If the queue is empty, this thread should block until another thread enqueues a value.
+ * This function should be concurrency-safe:
+ * multiple threads may call queue_enqueue() or queue_dequeue() simultaneously.
+ *
+ * @param queue the queue to remove from
+ * @return the value at the front of the queue
+ */
+void *queue_dequeue(queue_t *queue);
+
+/**
+ * Frees all resources associated with a heap-allocated queue.
+ * You may assume that the queue is already empty.
+ *
+ * @param queue a queue returned from queue_init()
+ */
+void queue_free(queue_t *queue);
+
+bool queue_ready(queue_t *queue);
+
+void queue_lock(queue_t *queue);
+
+void queue_unlock(queue_t *queue);
+
+void queue_sleep(queue_t *queue, bool condition);
+
+void queue_wake(queue_t *queue, bool condition);
+
+#endif /* QUEUE_H */

thread_pool.c

+#include "thread_pool.h"
+
+typedef struct thread_pool {
+    queue_t *job_queue;
+    pthread_t *workers;
+    size_t size;
+    bool finished;
+} thread_pool_t;
+
+typedef struct job {
+    work_function_t function;
+    void *aux;
+} job_t;
+
+void *get_do_work(void *aux) {
+    thread_pool_t *pool = (thread_pool_t *) aux;
+    queue_t *jobs = pool->job_queue;
+    while (true) {
+        job_t *new_job = queue_dequeue(jobs);
+        if (new_job->function == NULL) {
+            free(new_job);
+            break;
+        }
+        (*new_job->function)(new_job->aux);
+        free(new_job);
+    }
+    return NULL;
+}
+
+bool pool_ready(thread_pool_t *pool) {
+    return queue_ready(pool->job_queue);
+}
+
+thread_pool_t *thread_pool_init(size_t num_worker_threads) {
+    thread_pool_t *thread_pool = malloc(sizeof(thread_pool_t));
+    thread_pool->job_queue = queue_init();
+    thread_pool->workers = malloc((num_worker_threads) * sizeof(pthread_t));
+    thread_pool->size = num_worker_threads;
+    thread_pool->finished = false;
+    for (size_t i = 0; i < num_worker_threads; i++) {
+        pthread_create(&thread_pool->workers[i], NULL, get_do_work, thread_pool);
+    }
+    return thread_pool;
+}
+
+void thread_pool_add_work(thread_pool_t *pool, work_function_t function, void *aux) {
+    job_t *job = malloc(sizeof(job_t));
+    *job = (job_t){function, aux};
+    queue_enqueue(pool->job_queue, job);
+}
+
+void thread_pool_finish(thread_pool_t *pool) {
+    size_t size = pool->size;
+    for (size_t i = 0; i < size; i++) {
+        thread_pool_add_work(pool, NULL, NULL);
+    }
+    for (size_t i = 0; i < size; i++) {
+        pthread_join(pool->workers[i], NULL);
+    }
+    free(pool->workers);
+    queue_free(pool->job_queue);
+    free(pool);
+}
\ No newline at end of file

thread_pool.h

+#ifndef THREAD_POOL_H
+#define THREAD_POOL_H
+
+#include <pthread.h>
+#include <stdbool.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+#include "queue.h"
+
+/**
+ * A pool of threads which perform work in parallel.
+ * The pool contains a fixed number of threads specified in thread_pool_init()
+ * and a shared queue of work for the worker threads to run.
+ * Each worker thread dequeues new work from the queue when its current work is finished.
+ */
+typedef struct thread_pool thread_pool_t;
+
+/** A function that can run on a thread in a thread pool */
+typedef void (*work_function_t)(void *aux);
+
+/**
+ * Creates a new heap-allocated thread pool with the given number of worker threads.
+ * All worker threads should start immediately so they can perform work
+ * as soon as thread_pool_add_work() is called.
+ *
+ * @param num_worker_threads the number of threads in the pool
+ * @return a pointer to the new thread pool
+ */
+thread_pool_t *thread_pool_init(size_t num_worker_threads);
+
+/**
+ * Adds work to a thread pool.
+ * The work will be performed by a worker thread as soon as all previous work is finished.
+ *
+ * @param pool the thread pool to perform the work
+ * @param function the function to call on a thread in the thread pool
+ * @param aux the argument to call the work function with
+ */
+void thread_pool_add_work(thread_pool_t *pool, work_function_t function, void *aux);
+
+/**
+ * Waits for all work added to a thread pool to finish,
+ * then frees all resources associated with a heap-allocated thread pool.
+ * A special value (e.g. NULL) can be put in the work queue to mark the end of the work.
+ * thread_pool_add_work() cannot be used on this pool once this function is called.
+ *
+ * @param pool the thread pool to close
+ */
+void thread_pool_finish(thread_pool_t *pool);
+
+#endif /* THREAD_POOL_H */

wolff_cluster.c

+#include <stdio.h>
+#include <stdint.h>
+#include <math.h>
+#include "thread_pool.h"
+#include "pthread.h"
+#include "npy.h"
+/**
+ * CHANGE THE TEST VARIABLES AS NEEDED.
+ *
+ */
+#define N_LATTICES 1000
+#define N 128
+#define EPOCHS 20000
+#define FILE_NAME "data/2D_Ising_N%d_epochs_%d_beta_%g_lattices_%d.npy"
+
+
+const double BETA_VALS[] = {.39, .3925, .395, .3975, .4, .4025, .415, .4175, .42, .4225, .425, .4275, .43, .4325, .435, .4375, .44};
+
+
+/**
+ * DO NOT MESS WITH ANYTHING BELOW
+ *
+ */
+#define DIM 3
+
+const int  SHAPE[3] = {N_LATTICES, N, N};
+const size_t NUM_THREADS = 16;
+
+typedef struct point{int x, y;} point_t;
+
+
+
+int  **initialize_lattice()
+{
+    int  **s = calloc(N, sizeof(int  *));
+    for (size_t i = 0; i < N; i++)
+    {
+        s[i] = calloc(N, sizeof(int));
+        for (size_t j = 0; j < N; j++)
+        {
+            s[i][j] = (float)rand() / RAND_MAX < 0.5 ? +1 : -1;
+        }
+    }
+    return s;
+}
+
+bool **initialize_cluster()
+{
+    bool **s = calloc(N, sizeof(bool *));
+    for (size_t i = 0; i < N; i++)
+    {
+        s[i] = calloc(N, sizeof(bool));
+        for (size_t j = 0; j < N; j++)
+        {
+            s[i][j] = false;
+        }
+    }
+    return s;
+}
+
+void free_lattice(int  **lattice)
+{
+    for (size_t i = 0; i < N; i++)
+    {
+        free(lattice[i]);
+    }
+    free(lattice);
+}
+
+void free_cluster(bool **lattice)
+{
+    for (size_t i = 0; i < N; i++)
+    {
+        free(lattice[i]);
+    }
+    free(lattice);
+}
+
+point_t random_point()
+{
+    return (point_t){.x = rand() % N, .y = rand() % N};
+}
+
+point_t *get_nearest_neighbors(point_t site_index)
+{
+    point_t *nearest_neighbors = calloc(4, sizeof(point_t));
+    nearest_neighbors[0] = (point_t){.x = site_index.x == 0 ? N - 1 : site_index.x - 1, .y = site_index.y};
+    nearest_neighbors[1] = (point_t){.x = site_index.x == N - 1 ? 0 : site_index.x + 1, .y = site_index.y};
+    nearest_neighbors[2] = (point_t){.x = site_index.x, .y = site_index.y == 0 ? N - 1 : site_index.y + 1};
+    nearest_neighbors[3] = (point_t){.x = site_index.x, .y = site_index.y == N - 1 ? 0 : site_index.y + 1};
+    return nearest_neighbors;
+}
+
+void grow_cluster(point_t p, int cluster_spin, bool **cluster, int **s, double prob)
+{
+    point_t *nn = get_nearest_neighbors(p);
+    cluster[p.x][p.y] = true;
+    s[p.x][p.y] = -s[p.x][p.y];
+    for (size_t i = 0; i < 4; i++)
+    {
+        if (!cluster[nn[i].x][nn[i].y] && s[nn[i].x][nn[i].y] == cluster_spin && (double)rand() / (double)RAND_MAX  < prob)
+        {
+            grow_cluster(nn[i], cluster_spin, cluster, s, prob);
+        }
+    }
+}
+
+void one_monte_carlo_step(int **lattice, double beta, double k1)
+{
+    bool **cluster = initialize_cluster();
+    double p = 1 - exp(-2 * beta * k1);
+    point_t root = random_point();
+    int  spin = lattice[root.x][root.y];
+    grow_cluster(root, spin, cluster, lattice, p);
+    free_cluster(cluster);
+}
+
+void lattices_to_numpy(int ***data, double beta, int n, int n_lattices, int epochs){
+
+    char file_name[100];
+    sprintf(file_name, FILE_NAME, n, epochs, beta, n_lattices);
+    FILE *f = fopen(file_name, "wb");
+    skip_npy_header(f);
+    for (size_t j = 0; j < n_lattices; j++)
+    {
+        for (size_t k = 0; k < n; k++)
+        {
+            fwrite(data[j][k], sizeof(int), n, f);
+        }
+    }
+    int shape[] = {n_lattices, n, n};
+    write_npy_header(f, DIM, shape, NPY_INT);
+    fclose(f);
+}
+
+// extern void wolff_cluster_mc(float_t beta, float_t k1, int n, int n_lattices, int epochs){
+//     int ***data = calloc(n_lattices, sizeof(int **));
+//     double total_time = 0.0;
+//     for (size_t i = 0; i < n_lattices; i++)
+//     {
+//         clock_t begin = clock();
+//         int **lattice = initialize_lattice();
+//         for(size_t j = 0; j < epochs; j++){
+//             one_monte_carlo_step(lattice, beta, k1);
+//         }
+//         data[i] = lattice;
+//         clock_t end = clock();
+//         double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
+//         total_time += time_spent;
+//         printf("Current Beta Value: %lf\tCurrent Lattice: %d\tTime Spent: %lf\tTotal Time Elapsed: %lf\n", beta, i + 1, time_spent, total_time);
+//     }
+//     lattices_to_numpy(data, beta, n, n_lattices, epochs);
+//     printf("Overall Time: %lf\t Beta Value: %lf\n", total_time, beta);
+// }
+
+static void wolff_cluster_pthread(void *arg){
+    double beta = *(double*)arg;
+    int ***data = calloc(N_LATTICES, sizeof(int **));
+    int  **lattice = initialize_lattice();
+    for(size_t j = 0; j < N_LATTICES; j++){
+        for(size_t i = 0; i < EPOCHS; i++){
+            one_monte_carlo_step(lattice, beta, 1.0);
+        }
+        data[j] = lattice;
+    }
+    lattices_to_numpy(data, beta, N, N_LATTICES, EPOCHS);
+}
+
+int main(int argc, char *argv[])
+{
+    // int beta_num = sizeof(BETA_VALS) / sizeof(double);
+    // double total_time = 0.0;
+    // for (size_t i = 0; i < beta_num; i++)
+    // {
+    //     int ***data = calloc(N_LATTICES, sizeof(int **));
+    //     for (size_t j = 0; j < N_LATTICES; j++)
+    //     {
+    //         clock_t begin = clock();
+    //         double beta = BETA_VALS[i];
+    //         int **lattice = initialize_lattice();
+    //         for (size_t k = 0; k < EPOCHS; k++)
+    //         {
+    //             one_monte_carlo_step(lattice, beta, 1.0);
+    //         }
+    //         data[j] = lattice;
+    //         clock_t end = clock();
+    //         double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
+    //         total_time += time_spent;
+    //         printf("Current Beta Value: %lf\tCurrent Lattice: %d\tTime Spent: %lf\tTotal Time Elapsed: %lf\n", BETA_VALS[i], j, time_spent, total_time);
+    //     }
+    //     lattices_to_numpy(data, BETA_VALS[i], N, N_LATTICES, EPOCHS);
+    // }
+    // printf("Overall Time: %lf", total_time);
+
+    int beta_num = sizeof(BETA_VALS) / sizeof(double);
+    thread_pool_t *thread_pool = thread_pool_init(NUM_THREADS);
+    for(size_t i = 0; i < beta_num; i++){
+        thread_pool_add_work(thread_pool, wolff_cluster_pthread, (void*)&BETA_VALS[i]);
+    }
+    thread_pool_finish(thread_pool);
+}

wolff_cluster.h

+#ifndef WOLFF_CLUSTER_H 
+#define WOLFF_CLUSTER_H
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <stdbool.h>
+#include <string.h>
+#include <assert.h>
+#include <stdint.h>
+#include <pthread.h>
+#include "npy.h"
+#include "list.h"
+
+
+// extern void wolff_cluster_mc(float_t beta, float_t k1, int n, int n_lattices, int epochs);
+
+#endif

wolff_cluster/wolff_cluster.c

@@ -19,17 +19,17 @@ const double BETA_VALS[] = {.39, .3925, .395, .3975, .4, .4025, .415, .4175, .42
  */
 #define DIM 3
 
-const int SHAPE[3] = {N_LATTICES, N, N};
+const int64_t SHAPE[3] = {N_LATTICES, N, N};
 
 typedef struct point{int x, y;} point_t;
 
 int **initialize_lattice()
 {
-    int **s = calloc(N, sizeof(int *));
-    for (int i = 0; i < N; i++)
+    int64_t **s = calloc(N, sizeof(int64_t *));
+    for (size_t i = 0; i < N; i++)
     {
         s[i] = calloc(N, sizeof(int));
-        for (int j = 0; j < N; j++)
+        for (size_t j = 0; j < N; j++)
         {
             s[i][j] = (float)rand() / RAND_MAX < 0.5 ? +1 : -1;
         }
@@ -40,10 +40,10 @@ int **initialize_lattice()
 bool **initialize_cluster()
 {
     bool **s = calloc(N, sizeof(bool *));
-    for (int i = 0; i < N; i++)
+    for (size_t i = 0; i < N; i++)
     {
         s[i] = calloc(N, sizeof(bool));
-        for (int j = 0; j < N; j++)
+        for (size_t j = 0; j < N; j++)
         {
             s[i][j] = false;
         }
@@ -51,9 +51,9 @@ bool **initialize_cluster()
     return s;
 }
 
-void free_lattice(int **lattice)
+void free_lattice(int64_t **lattice)
 {
-    for (int i = 0; i < N; i++)
+    for (size_t i = 0; i < N; i++)
     {
         free(lattice[i]);
     }
@@ -62,7 +62,7 @@ void free_lattice(int **lattice)
 
 void free_cluster(bool **lattice)
 {
-    for (int i = 0; i < N; i++)
+    for (size_t i = 0; i < N; i++)
     {
         free(lattice[i]);
     }
@@ -89,7 +89,7 @@ void grow_cluster(point_t p, int cluster_spin, bool **cluster, int **s, double p
     point_t *nn = get_nearest_neighbors(p);
     cluster[p.x][p.y] = true;
     s[p.x][p.y] = -s[p.x][p.y];
-    for (int i = 0; i < 4; i++)
+    for (size_t i = 0; i < 4; i++)
     {
         if (!cluster[nn[i].x][nn[i].y] && s[nn[i].x][nn[i].y] == cluster_spin && (double)rand() / (double)RAND_MAX  < prob)
         {
@@ -103,7 +103,7 @@ void one_monte_carlo_step(int **lattice, double beta, double k1)
     bool **cluster = initialize_cluster();
     double p = 1 - exp(-2 * beta * k1);
     point_t root = random_point();
-    int spin = lattice[root.x][root.y];
+    int64_t spin = lattice[root.x][root.y];
     grow_cluster(root, spin, cluster, lattice, p);
     free_cluster(cluster);
 }
@@ -114,9 +114,9 @@ void lattices_to_numpy(int ***data, float_t beta, int n, int n_lattices, int epo
     sprintf(file_name, FILE_NAME, n, epochs, beta, n_lattices);
     FILE *f = fopen(file_name, "wb");
     skip_npy_header(f);
-    for (int j = 0; j < n_lattices; j++)
+    for (size_t j = 0; j < n_lattices; j++)
     {
-        for (int k = 0; k < n; k++)
+        for (size_t k = 0; k < n; k++)
         {
             fwrite(data[j][k], sizeof(int), n, f);
         }
@@ -129,11 +129,11 @@ void lattices_to_numpy(int ***data, float_t beta, int n, int n_lattices, int epo
 extern void wolff_cluster_mc(float_t beta, float_t k1, int n, int n_lattices, int epochs){
     int ***data = calloc(n_lattices, sizeof(int **));
     double total_time = 0.0;
-    for (int i = 0; i < n_lattices; i++)
+    for (size_t i = 0; i < n_lattices; i++)
     {
         clock_t begin = clock();
         int **lattice = initialize_lattice();
-        for(int j = 0; j < epochs; j++){
+        for(size_t j = 0; j < epochs; j++){
             one_monte_carlo_step(lattice, beta, k1);
         }
         data[i] = lattice;
@@ -150,15 +150,15 @@ int main(int argc, char *argv[])
 {
     int beta_num = sizeof(BETA_VALS) / sizeof(double);
     double total_time = 0.0;
-    for (int i = 0; i < beta_num; i++)
+    for (size_t i = 0; i < beta_num; i++)
     {
         int ***data = calloc(N_LATTICES, sizeof(int **));
-        for (int j = 0; j < N_LATTICES; j++)
+        for (size_t j = 0; j < N_LATTICES; j++)
         {
             clock_t begin = clock();
             double beta = BETA_VALS[i];
             int **lattice = initialize_lattice();
-            for (int k = 0; k < EPOCHS; k++)
+            for (size_t k = 0; k < EPOCHS; k++)
             {
                 one_monte_carlo_step(lattice, beta, 1.0);
             }

wolff_cluster/wolff_cluster.h

@@ -8,6 +8,7 @@
 #include <stdbool.h>
 #include <string.h>
 #include <assert.h>
+#include <stdint.h>
 #include "npy.h"