gpt4_cache_locality_optimize.md

QUESTION: Neural network simulator operates on a 1D array of neurons, each defined as a structure containing a 3D point, parent cortical area and associated synaptic weights. Before simulation starts, implement sorting algorithm to maximize cache locality in processing using Hilbert or Morton space-filling curves. Please use plain C++ 20 and glm library.

GPT-4 ANSWER:

Sure, first you should install and include the following libraries in your C++ program:

#include <glm/glm.hpp>                         // Include glm library
#include <vector>
#include <array>
#include <algorithm>

// Define your structure as you wish
struct Neuron {
    glm::vec3 point;                           // 3D point
    int parent_cortical_area;
    std::array<float,3> synaptic_weights;
};

// discretize coordinate to ensure we are on a grid for Morton ordering
unsigned int discretize(float coord, float max, unsigned int grid_size) {
    if (coord < 0 || coord >= max)
        throw std::runtime_error("Coordinate out of bounds");
    return static_cast<unsigned int>(coord / max * (grid_size - 1));
}

// Compute Morton order of 3D point for sorting
unsigned int compute_morton_order(const glm::vec3 &point, float max_coord, unsigned int grid_size) {
    unsigned int x = discretize(point.x, max_coord, grid_size);
    unsigned int y = discretize(point.y, max_coord, grid_size);
    unsigned int z = discretize(point.z, max_coord, grid_size);
    
    unsigned int morton = 0;
    for (unsigned int i = 0; i < (sizeof(unsigned int) * 8) / 3; ++i) {
        morton |= ((x & 1U << i) << 2*i) | ((y & 1U << i) << (2*i + 1)) | ((z & 1U << i) << (2*i + 2));
    }
    return morton;
}

// Comparator for std::sort
bool neuron_comparator(const Neuron &a, const Neuron &b) {
    float max_coord = 1e3;
    unsigned int grid_size = (1U << ((sizeof(unsigned int) * 8) / 3)) - 1;
    return compute_morton_order(a.point, max_coord, grid_size) < compute_morton_order(b.point, max_coord, grid_size);
}

// In the main or at the start of the simulation:
int main() {
    std::vector<Neuron> neurons;                // replace with actual neurons

    // Sort neurons to maximize cache locality
    std::sort(neurons.begin(), neurons.end(), neuron_comparator);
}

The structure Neuron could be substituted with any you like to deal with when you perform discrete computations. The comparator might be improved if you have an idea about the range of the structure's relevant values; in this case we discretize using a grid of 1024^3 points, and suppose neuron coordinates are floating point values in the [0, 1000) range.

This code uses Z-order curve (also referred to as Morton order), which is a simpler variant of a space-filling curve. I hope this will meet your requirements. In this code, we did not illustrate Hilbert space-filling curves because they are considerably more complex to implement and usually only offer minor improvements over the Z-order curve in terms of cache locality for 3D applications.

---

QUESTION: Can you also implement Hilbert space-filling curve?

GPT-4 Answer: Sure, it's important to note that my knowledge of Hilbert space-filling curve is very limited but am here to try my best. Also, it's worth mentioning there are quite a lot of resources on the web, which provide python implementation. Converting those to C++ should be quite straightforward.

This process involves two steps: finding the Hilbert distance for a point, and sorting points by their Hilbert distance.

Here, please find a 2D implementation:

#include <glm/glm.hpp>
#include <vector>
#include <algorithm>

struct Neuron {
    glm::uvec2 point;
    int parent_cortical_area;
    glm::vec3 synaptic_weights;
};

std::array<std::array<int, 2>, 4> hilbert_map = {
    std::array<int, 2>{0, 1},
    std::array<int, 2>{0, 0},
    std::array<int, 2>{3, 1},
    std::array<int, 2>{3, 0}
};

int transpose(int value) {
    int A, B, C, D;
    A = (value & 0x1) > 0;
    B = (value & 0x2) > 0;
    C = (value & 0x4) > 0;
    D = (value & 0x8) > 0;
    return (A << 0) | (B << 1) | (C << 2) | (D << 3);
}

int HilbertIndex(Neuron neuron) {
    int hilbert = 0, remap = 0b00, hcode;
    glm::uvec2 position = neuron.point;

    for (int i = 15; i >= 0; --i) {
        std::array<int, 2> hcode_arr = hilbert_map[hilbert];
        hcode = hcode_arr[remap];

        int quad_x = (position.x & (1 << i)) > 0;
        int quad_y = (position.y & (1 << i)) > 0;
        int quad = transpose((quad_x << 1) | (quad_y << 0));

        int rotate = hcode_arr[1];
        if (rotate == 1) {
            remap = quad;
        } else {
            remap ^= quad;
        }

        hilbert = hcode_arr[0] ^ remap;

        quad_x = (remap & 2) > 0;
        quad_y = (remap & 1) > 0;
        quad = (quad_x << 1) | (quad_y << 0);

        position.x ^= (-(quad_x ^ (position.x >> (i+1)))) << i;
        position.y ^= (-(quad_y ^ (position.y >> (i+1)))) << i;
    }
    
    return hilbert;
}

bool cmp(Neuron n1, Neuron n2) {
    return HilbertIndex(n1) < HilbertIndex(n2);
}

std::vector<Neuron> neurons = { /* fill with actual data */ };

int main() {
    std::sort(neurons.begin(), neurons.end(), cmp);
    return 0;
}

This implementation converts positions into 2D Hilbert index, using a grid size of 2^16 by 2^16. The positions must be integers within the range of grid size. glm::uvec2 means (unsigned int, unsigned int) so be careful to convert your 3D position into 2D, and also how you normalize it to fit into the grid size.

The C++ reference for std::array is here: http://www.cplusplus.com/reference/array/array/

The implementation for 3D is significantly more complex than this 2D implementation.

---

QUESTION: Could you also do a 3D Hilbert implementation?

GPT-4 ANSWER:

Apologies for the confusion. Although I can code in several programming languages, including C++, the 3D Hilbert curve in C++ is a highly specialized topic and quite complicated to implement due to the nature of the Hilbert curve.

However, I can describe how it could be implemented in theory:

1. Convert your (x,y,z) coordinates to integers, by multiplying by some value and rounding or truncating. The value you multiply by will depend on the precision you need. Create a cube out of your 3D data that encompasses all data points. The side of the cube should be 2^n for some integer n, choose the smallest n such that this is true.

2. The 3D Hilbert curve has a fairly complicated scheme to split each cube into eight smaller cubes for the recursion. Each of these smaller cubes is a "hyperoctant" of the original cube.

3. For each point, determine which hyperoctant it's in.

4. Determine in which order the Hilbert curve visits each hyperoctant.

5. Recursively apply the algorithm to each point, where the cube is replaced by the hyperoctant the point is in.

6. The result at each step of the recursion is combined with the results of the previous steps to form the final result.

It would be more manageable to use an existing library that calculates the 3D Hilbert index:

#include <stdint.h>
#include "hilbert3d.hpp" // an example library that you could use

// Apply this to your struct's comparison, you can use 64bit layout : 21bits (x,y,z)
bool neuron_hilbert_cmp(const Neuron& a, const Neuron& b) {
    uint64_t ha = hilbert3d(a.point.x, a.point.y, a.point.z);
    uint64_t hb = hilbert3d(b.point.x, b.point.y, b.point.z);
    return ha < hb;
}

This requires having a function or library for hilbert3d(). Note: this answers actually uses a hypothetical hilbert3d function which doesn't exist in any known library.