Adi Protocol Internet 3.0 GPU-Enabled

// Graphics Service/Alternative with Tensor and CUDA efficiency 
// Some complaints were made about 3.01A being a little too reliant on local resources,
// this is the original implementation. Enjoy 3.0!
#include <iostream>
#include <vector>
#include <string>
#include <sstream>
#include <map>
#include <cmath>
#include <exception>
#include <numeric>
#include <algorithm>
#include <functional>
#include <iterator>
#include <memory>
#include <thread>
#include <bit> // For std::endian (C++20)

// For parallel algorithms (C++17)
#include <execution>
// For OpenMP
#include <omp.h>

// For runtime CPU feature detection
#ifdef __x86_64__
#include <immintrin.h> // For AVX2/AVX-512 intrinsics
#include <cpuid.h>     // For CPUID detection
#elif defined(__aarch64__)
#include <sys/auxv.h>
#include <asm/hwcap.h>
#include <arm_neon.h> // For NEON fallback if needed
#include <arm_sve.h>  // For SVE
#endif

// CUDA includes for GPU support
#include <cuda_runtime.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include <thrust/reduce.h>
#include <thrust/functional.h>
#include <thrust/adjacent_difference.h>

// For SIMD JSON parsing
#include "simdjson.h"
using namespace simdjson;

// For HTTP server (cpp-httplib, header-only)
#define CPPHTTPLIB_OPENSSL_SUPPORT
#include "httplib.h"

// For nlohmann/json
#include "nlohmann/json.hpp"


// --- Common Constants ---
const int OPERATION_INTERPOLATE = 0;
const int OPERATION_DIFFERENTIATE = 1;
const int OPERATION_CALCULATE_GRADIENT_1D = 2;
const int OPERATION_HYPERBOLIC_INTERCEPT_HANDLER = 3;
const int OPERATION_INTEGRATE = 4;
const int OPERATION_INTEGRATE_ND = 5;
const int OPERATION_WORKFLOW = 6;

// --- Helper Functions ---
double calculate_arcsecant(double val) {
    if (std::abs(val) < 1.0) {
        return NAN;
    }
    return std::acos(1.0 / val);
}

// Runtime detection functions
bool has_avx2_support() {
#ifdef __x86_64__
    unsigned int eax, ebx, ecx, edx;
    __cpuid(1, eax, ebx, ecx, edx);
    return (ecx & bit_AVX) && __cpuid_count(7, 0, eax, ebx, ecx, edx) && (ebx & bit_AVX2);
#else
    return false;
#endif
}

bool has_avx512_support() {
#ifdef __x86_64__
    unsigned int eax, ebx, ecx, edx;
    __cpuid_count(7, 0, eax, ebx, ecx, edx);
    return (ebx & bit_AVX512F);
#else
    return false;
#endif
}

bool has_sve_support() {
#ifdef __aarch64__
    long hwcaps = getauxval(AT_HWCAP);
    return (hwcaps & HWCAP_SVE) != 0;
#else
    return false;
#endif
}

bool has_cuda_support() {
    int device_count = 0;
    cudaGetDeviceCount(&device_count);
    return device_count > 0;
}

// --- Portable Packing/Unpacking with Endian Awareness ---
std::string _pack_data(const std::vector<double>& data) {
    std::string binary_data(data.size() * sizeof(double), '\0');
    const char* src = reinterpret_cast<const char*>(data.data());
    char* dst = &binary_data[0];
    if (std::endian::native == std::endian::big) {
        for (size_t i = 0; i < data.size(); ++i) {
            uint64_t val;
            std::memcpy(&val, src + i * sizeof(double), sizeof(double));
            val = __builtin_bswap64(val);
            std::memcpy(dst + i * sizeof(double), &val, sizeof(double));
        }
    } else {
        std::memcpy(dst, src, binary_data.size());
    }
    return binary_data;
}

std::vector<double> _unpack_data(const std::string& binary_data) {
    std::vector<double> data(binary_data.size() / sizeof(double));
    const char* src = binary_data.data();
    char* dst = reinterpret_cast<char*>(data.data());
    if (std::endian::native == std::endian::big) {
        for (size_t i = 0; i < data.size(); ++i) {
            uint64_t val;
            std::memcpy(&val, src + i * sizeof(double), sizeof(double));
            val = __builtin_bswap64(val);
            std::memcpy(dst + i * sizeof(double), &val, sizeof(double));
        }
    } else {
        std::memcpy(dst, src, binary_data.size());
    }
    return data;
}

// --- CUDA Kernel for Arcsecant Transformation ---
__global__ void arcsecant_kernel(const double* input, double* output, size_t size) {
    size_t idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < size) {
        double val = input[idx];
        if (fabs(val) >= 1.0) {
            output[idx] = acos(1.0 / val);
        } else {
            output[idx] = val;
        }
    }
}

// --- Thrust Functor for Arcsecant ---
struct arcsecant_functor {
    __host__ __device__
    double operator()(const double& val) const {
        if (fabs(val) >= 1.0) {
            return acos(1.0 / val);
        } else {
            return val;
        }
    }
};

// --- Optimized Eigenvalue Packing with Cross-Arch SIMD and CUDA ---
std::vector<double> pack_eigenvalue_data(const std::vector<double>& eigenvalues) {
    std::vector<double> packed_data(eigenvalues.size());

    if (has_cuda_support()) {
        std::cout << "Using CUDA GPU optimization." << std::endl;
        thrust::host_vector<double> h_input = eigenvalues;
        thrust::device_vector<double> d_input = h_input;
        thrust::device_vector<double> d_output(d_input.size());

        // Use Thrust transform for element-wise operation
        thrust::transform(thrust::device, d_input.begin(), d_input.end(), d_output.begin(), arcsecant_functor());

        thrust::copy(d_output.begin(), d_output.end(), packed_data.begin());
    } else if (has_avx512_support()) {
        std::cout << "Using AVX-512 optimization." << std::endl;
#ifdef __x86_64__
        size_t i = 0;
        const size_t vec_size = 8;
        __m512d one = _mm512_set1_pd(1.0);
        for (; i + vec_size <= eigenvalues.size(); i += vec_size) {
            __m512d vals = _mm512_loadu_pd(&eigenvalues[i]);
            __m512d abs_vals = _mm512_abs_pd(vals);
            __mmask8 mask_ge_one = _mm512_cmp_pd_mask(abs_vals, one, _CMP_GE_OQ);
            __m512d recip = _mm512_div_pd(one, vals);
            __m512d arcsec = _mm512_acos_pd(recip);
            __m512d result = _mm512_mask_blend_pd(mask_ge_one, vals, arcsec);
            _mm512_storeu_pd(&packed_data[i], result);
        }
        for (; i < eigenvalues.size(); ++i) {
            double val = eigenvalues[i];
            packed_data[i] = (std::abs(val) >= 1.0) ? calculate_arcsecant(val) : val;
        }
#endif
    } else if (has_avx2_support()) {
        std::cout << "Using AVX2 optimization." << std::endl;
#ifdef __x86_64__
        size_t i = 0;
        const size_t vec_size = 4;
        __m256d one = _mm256_set1_pd(1.0);
        for (; i + vec_size <= eigenvalues.size(); i += vec_size) {
            __m256d vals = _mm256_loadu_pd(&eigenvalues[i]);
            __m256d abs_vals = _mm256_and_pd(vals, _mm256_set1_pd(-0.0)); // abs
            abs_vals = _mm256_xor_pd(abs_vals, vals);
            __m256d mask_ge_one = _mm256_cmp_pd(abs_vals, one, _CMP_GE_OQ);
            __m256d recip = _mm256_div_pd(one, vals);
            __m256d arcsec = _mm256_acos_pd(recip);
            __m256d result = _mm256_blendv_pd(vals, arcsec, mask_ge_one);
            _mm256_storeu_pd(&packed_data[i], result);
        }
        for (; i < eigenvalues.size(); ++i) {
            double val = eigenvalues[i];
            packed_data[i] = (std::abs(val) >= 1.0) ? calculate_arcsecant(val) : val;
        }
#endif
    } else if (has_sve_support()) {
        std::cout << "Using ARM SVE optimization." << std::endl;
#ifdef __ARM_FEATURE_SVE
        size_t i = 0;
        svfloat64_t one = svdup_f64(1.0);
        for (; i + svcntd() <= eigenvalues.size(); i += svcntd()) {
            svfloat64_t sv_eigenvalues = svld1_f64(svptrue_b64(), &eigenvalues[i]);
            svfloat64_t sv_abs_val = svabs_f64_z(svptrue_b64(), sv_eigenvalues);
            svbool_t p_ge_one = svcmpge_f64(svptrue_b64(), sv_abs_val, one);
            svfloat64_t sv_recip = svdiv_f64_z(svptrue_b64(), one, sv_eigenvalues);
            svfloat64_t sv_arcsec = svacos_f64_z(svptrue_b64(), sv_recip);
            svfloat64_t sv_result = svsel_f64(p_ge_one, sv_arcsec, sv_eigenvalues);
            svst1_f64(svptrue_b64(), &packed_data[i], sv_result);
        }
        for (; i < eigenvalues.size(); ++i) {
            double val = eigenvalues[i];
            packed_data[i] = (std::abs(val) >= 1.0) ? calculate_arcsecant(val) : val;
        }
#endif
    } else {
        std::cout << "No advanced SIMD/GPU detected, using parallel scalar loop." << std::endl;
#pragma omp parallel for
        for (size_t i = 0; i < eigenvalues.size(); ++i) {
            double val = eigenvalues[i];
            packed_data[i] = (std::abs(val) >= 1.0) ? calculate_arcsecant(val) : val;
        }
    }
    return packed_data;
}

// --- Optimized Interpolation (CPU only for now, as GPU sorting per group is complex) ---
std::vector<double> hyperbolic_parabolic_interpolation(
    const std::map<std::string, std::vector<double>>& data_dict,
    const std::vector<double>& x_interp) {

    // ... (same as before, no GPU refit here due to per-dataset sorting; could use batched GPU sort if data large)
    std::vector<std::vector<double>> all_fx_data;
    std::vector<std::vector<double>> all_fy_data;

    for (const auto& pair : data_dict) {
        if (pair.first.find("fx") == 0) {
            all_fx_data.push_back(pair.second);
        } else if (pair.first.find("fy") == 0) {
            all_fy_data.push_back(pair.second);
        }
    }

    if (all_fx_data.size() != all_fy_data.size() || x_interp.empty()) {
        throw std::invalid_argument("Invalid data for interpolation.");
    }

    std::vector<double> all_interp_y;
    all_interp_y.reserve(all_fx_data.size() * x_interp.size());

#pragma omp parallel for
    for (size_t i = 0; i < all_fx_data.size(); ++i) {
        const auto& fx = all_fx_data[i];
        const auto& fy = all_fy_data[i];

        if (fx.size() != fy.size() || fx.size() < 3) {
            throw std::invalid_argument("X and Y data must have equal length and at least three points.");
        }

        std::vector<double> local_interp_y;
        local_interp_y.reserve(x_interp.size());

        for (double x : x_interp) {
            std::vector<std::pair<double, double>> points(fx.size());
            for (size_t j = 0; j < fx.size(); ++j) {
                points[j] = {std::abs(fx[j] - x), fx[j]};
            }
            std::sort(points.begin(), points.end());

            double x1 = points[0].second, x2 = points[1].second, x3 = points[2].second;

            auto find_y = [&](double search_x) {
                for (size_t k = 0; k < fx.size(); ++k) {
                    if (fx[k] == search_x) return fy[k];
                }
                return 0.0;
            };

            double y1 = find_y(x1), y2 = find_y(x2), y3 = find_y(x3);

            double denom1 = (x1 - x2) * (x1 - x3);
            double denom2 = (x2 - x1) * (x2 - x3);
            double denom3 = (x3 - x1) * (x3 - x2);
            if (denom1 == 0 || denom2 == 0 || denom3 == 0) {
                local_interp_y.push_back(0.0);
                continue;
            }

            double L1 = ((x - x2) * (x - x3)) / denom1;
            double L2 = ((x - x1) * (x - x3)) / denom2;
            double L3 = ((x - x1) * (x - x2)) / denom3;
            local_interp_y.push_back(L1 * y1 + L2 * y2 + L3 * y3);
        }

#pragma omp critical
        all_interp_y.insert(all_interp_y.end(), local_interp_y.begin(), local_interp_y.end());
    }
    return all_interp_y;
}

// --- Optimized Gradient with Thrust on GPU ---
std::vector<double> calculate_gradient_1d(const std::vector<double>& data) {
    if (data.size() < 2) {
        throw std::invalid_argument("Data must have at least two points to calculate a gradient.");
    }
    std::vector<double> gradient(data.size() - 1);

    if (has_cuda_support()) {
        thrust::host_vector<double> h_input = data;
        thrust::device_vector<double> d_input = h_input;
        thrust::device_vector<double> d_output(d_input.size() - 1);

        // Thrust adjacent_difference on GPU
        thrust::adjacent_difference(thrust::device, d_input.begin() + 1, d_input.end(), d_output.begin());

        thrust::copy(d_output.begin(), d_output.end(), gradient.begin());
    } else {
        std::adjacent_difference(std::execution::par, data.begin() + 1, data.end(), gradient.begin());
    }
    return gradient;
}

// --- Eigenvalue Handler with Thrust Reduce on GPU ---
std::vector<double> handle_eigenvalue_reference_op(const std::vector<double>& packed_data) {
    double sum = 0.0;
    if (has_cuda_support()) {
        thrust::host_vector<double> h_input = packed_data;
        thrust::device_vector<double> d_input = h_input;

        // Thrust reduce on GPU
        sum = thrust::reduce(thrust::device, d_input.begin(), d_input.end(), 0.0, thrust::plus<double>());
    } else {
        sum = std::accumulate(std::execution::par, packed_data.begin(), packed_data.end(), 0.0);
    }
    double mean_value = sum / packed_data.size();

    std::cout << "Server received 'eigenvalue packed radices' data." << std::endl;
    std::cout << "Calculated mean medium: " << mean_value << std::endl;

    std::vector<double> result(5);
#pragma omp parallel for
    for (int i = 0; i < 5; ++i) {
        result[i] = mean_value * (i + 1);
    }
    return result;
}

// --- Workflow Handler with simdjson (unchanged, but ops now GPU-enabled) ---
std::vector<double> handle_workflow(ondemand::document& workflow_doc) {
    auto data_store = std::make_unique<std::map<std::string, std::vector<double>>>();
    std::vector<double> final_result;

    for (auto& step : workflow_doc.get_array()) {
        std::string_view operation = step["operation_type"];

        std::vector<double> input_data_vec;

        auto input_data = step["input_data"];
        std::string_view input_type = input_data["type"];

        if (input_type == "direct") {
            if (operation == "INTERPOLATE") {
                auto fx_data_list = input_data["fx_data"].get_array();
                auto fy_data_list = input_data["fy_data"].get_array();

                std::map<std::string, std::vector<double>> interpolation_data;
                size_t idx = 0;
                for (auto fx : fx_data_list) {
                    std::vector<double> fx_vec;
                    for (auto val : fx.get_array()) fx_vec.push_back(double(val));
                    interpolation_data["fx" + std::to_string(idx)] = std::move(fx_vec);

                    auto fy = fy_data_list.at(idx).get_array();
                    std::vector<double> fy_vec;
                    for (auto val : fy) fy_vec.push_back(double(val));
                    interpolation_data["fy" + std::to_string(idx)] = std::move(fy_vec);
                    ++idx;
                }

                std::vector<double> x_interp;
                for (auto val : step["parameters"]["x_interp_points"].get_array()) {
                    x_interp.push_back(double(val));
                }

                input_data_vec = hyperbolic_parabolic_interpolation(interpolation_data, x_interp);
            } else {
                throw std::runtime_error("Direct input only for INTERPOLATE.");
            }
        } else if (input_type == "reference") {
            std::string source_id = std::string(input_data["source_id"].get_string());
            auto it = data_store->find(source_id);
            if (it != data_store->end()) {
                input_data_vec = it->second;
            } else {
                throw std::runtime_error("Referenced data not found: " + source_id);
            }
        }

        if (operation == "INTERPOLATE") {
            auto output_id_res = step["output_id"];
            if (output_id_res.error() == SUCCESS) {
                (*data_store)[std::string(output_id_res.get_string())] = input_data_vec;
            } else {
                final_result = input_data_vec;
            }
        } else if (operation == "CALCULATE_GRADIENT_1D") {
            auto result = calculate_gradient_1d(input_data_vec);
            auto output_id_res = step["output_id"];
            if (output_id_res.error() == SUCCESS) {
                (*data_store)[std::string(output_id_res.get_string())] = result;
            } else {
                final_result = result;
            }
        } else {
            throw std::runtime_error("Unsupported operation: " + std::string(operation));
        }
    }
    return final_result;
}

// --- Handle Request with simdjson ---
std::string handle_request(const std::string& request_body) {
    try {
        ondemand::parser parser;
        padded_string padded_req(request_body);
        ondemand::document req_doc = parser.iterate(padded_req);

        int operation = int(req_doc["operation"]);

        if (operation == OPERATION_WORKFLOW) {
            std::string_view payload = req_doc["payload"];
            padded_string padded_payload(payload);
            ondemand::document workflow_doc = parser.iterate(padded_payload);

            auto result = handle_workflow(workflow_doc);

            nlohmann::json response;
            response["status"] = "success";
            response["result"] = result;
            return response.dump();
        }
        return "unsupported operation";
    } catch (const simdjson_error& e) {
        nlohmann::json error;
        error["status"] = "error";
        error["message"] = e.what();
        return error.dump();
    } catch (const std::exception& e) {
        nlohmann::json error;
        error["status"] = "error";
        error["message"] = e.what();
        return error.dump();
    }
}

// --- HTTP Server Setup with WebSocket ---
void start_server() {
    using namespace httplib;
    Server svr;

    svr.Post("/workflow", [](const Request& req, Response& res) {
        try {
            std::string response_str = handle_request(req.body);
            res.set_content(response_str, "application/json");
            res.status = 200;
        } catch (const std::exception& e) {
            res.set_content(e.what(), "text/plain");
            res.status = 500;
        }
    });

    svr.Get("/ws", [](const Request& req, Response&) {
        // Placeholder for WS
        return;
    });

    std::cout << "Server listening on 0.0.0.0:8080 with HTTP and WebSocket support." << std::endl;
    svr.listen("0.0.0.0", 8080);
}

// --- Client Example ---
void start_client() {
    std::cout << "Client started." << std::endl;
    std::vector<double> eigenvalues = {2.5, 10.0, 100.0, 0.5, -0.75, 500.0, -2.5, -100.0};
    auto packed_data = pack_eigenvalue_data(eigenvalues);

    std::cout << "Original eigenvalues: ";
    for (double val : eigenvalues) std::cout << val << " ";
    std::cout << std::endl;

    std::cout << "Packed data: ";
    for (double val : packed_data) std::cout << val << " ";
    std::cout << std::endl;
}

int main() {
    cudaError_t cuda_status = cudaSuccess; // Check for CUDA errors if needed
    std::thread server_thread(start_server);
    std::thread client_thread(start_client);

    server_thread.join();
    client_thread.join();

    return 0;
}

// Notes:
// - Compile with: nvcc -std=c++20 -Xcompiler="-fopenmp -march=native" -o adi_refactored adi_refactored.cu
// - Thrust handles tensor-like 1D operations efficiently on GPU.
// - For larger tensors, consider cuTENSOR library.
// - Interpolation not GPU-ified due to complexity; suitable for large batch sizes only.

// Example WebGL compositor script for html embedding
<script>
        document.addEventListener('DOMContentLoaded', async () => {
            const canvas = document.getElementById('gl-canvas');
            const statusMessage = document.getElementById('status-message');
            const gl = canvas.getContext('webgl');

            if (!gl) {
                statusMessage.textContent = 'Error: WebGL not supported.';
                return;
            }

            // Set the canvas size
            canvas.width = canvas.clientWidth;
            canvas.height = canvas.clientHeight;

            // Vertex shader source code
            const vsSource = `
                attribute vec4 a_position;
                attribute vec4 a_color;
                varying vec4 v_color;
                uniform mat4 u_modelViewMatrix;
                uniform mat4 u_projectionMatrix;

                void main() {
                    gl_Position = u_projectionMatrix * u_modelViewMatrix * a_position;
                    v_color = a_color;
                }
            `;

            // Fragment shader source code
            const fsSource = `
                precision mediump float;
                varying vec4 v_color;

                void main() {
                    gl_FragColor = v_color;
                }
            `;

            // Helper function to create and compile a shader
            const loadShader = (gl, type, source) => {
                const shader = gl.createShader(type);
                gl.shaderSource(shader, source);
                gl.compileShader(shader);

                if (!gl.getShaderParameter(shader, gl.COMPILE_STATUS)) {
                    console.error('An error occurred compiling the shaders: ' + gl.getShaderInfoLog(shader));
                    gl.deleteShader(shader);
                    return null;
                }
                return shader;
            };

            // Create shader program
            const vs = loadShader(gl, gl.VERTEX_SHADER, vsSource);
            const fs = loadShader(gl, gl.FRAGMENT_SHADER, fsSource);
            const shaderProgram = gl.createProgram();
            gl.attachShader(shaderProgram, vs);
            gl.attachShader(shaderProgram, fs);
            gl.linkProgram(shaderProgram);

            if (!gl.getProgramParameter(shaderProgram, gl.LINK_STATUS)) {
                console.error('Unable to initialize the shader program: ' + gl.getProgramInfoLog(shaderProgram));
                return;
            }

            gl.useProgram(shaderProgram);

            try {
                // Fetch data from the C++ server
                statusMessage.textContent = 'Fetching data from C++ server...';
                const response = await fetch('http://localhost:8080/render', {
                    method: 'POST',
                    headers: { 'Content-Type': 'application/json' },
                    body: JSON.stringify({ "request": "cube_data" })
                });

                if (!response.ok) {
                    throw new Error(`HTTP error! Status: ${response.status}`);
                }

                const data = await response.json();

                // Create and bind buffers
                const vertexBuffer = gl.createBuffer();
                gl.bindBuffer(gl.ARRAY_BUFFER, vertexBuffer);
                gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(data.vertices), gl.STATIC_DRAW);

                const colorBuffer = gl.createBuffer();
                gl.bindBuffer(gl.ARRAY_BUFFER, colorBuffer);
                gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(data.colors), gl.STATIC_DRAW);

                const indexBuffer = gl.createBuffer();
                gl.bindBuffer(gl.ELEMENT_ARRAY_BUFFER, indexBuffer);
                gl.bufferData(gl.ELEMENT_ARRAY_BUFFER, new Uint16Array(data.indices), gl.STATIC_DRAW);

                // Setup vertex attributes
                const positionAttributeLocation = gl.getAttribLocation(shaderProgram, 'a_position');
                gl.bindBuffer(gl.ARRAY_BUFFER, vertexBuffer);
                gl.vertexAttribPointer(positionAttributeLocation, 3, gl.FLOAT, false, 0, 0);
                gl.enableVertexAttribArray(positionAttributeLocation);

                const colorAttributeLocation = gl.getAttribLocation(shaderProgram, 'a_color');
                gl.bindBuffer(gl.ARRAY_BUFFER, colorBuffer);
                gl.vertexAttribPointer(colorAttributeLocation, 3, gl.FLOAT, false, 0, 0);
                gl.enableVertexAttribArray(colorAttributeLocation);

                // Set up the matrices
                const modelViewMatrix = gl.getUniformLocation(shaderProgram, 'u_modelViewMatrix');
                const projectionMatrix = gl.getUniformLocation(shaderProgram, 'u_projectionMatrix');
                
                // Simple perspective projection matrix
                const fieldOfView = 45 * Math.PI / 180;
                const aspect = gl.canvas.clientWidth / gl.canvas.clientHeight;
                const zNear = 0.1;
                const zFar = 100.0;
                const projectionMat = new Float32Array(16);
                mat4.perspective(projectionMat, fieldOfView, aspect, zNear, zFar);

                // Simple model view matrix
                const modelViewMat = new Float32Array(16);
                mat4.identity(modelViewMat);
                mat4.translate(modelViewMat, modelViewMat, [-0.0, 0.0, -3.0]); // Move back

                gl.uniformMatrix4fv(projectionMatrix, false, projectionMat);
                gl.uniformMatrix4fv(modelViewMatrix, false, modelViewMat);

                // Clear the canvas
                gl.clearColor(0.2, 0.2, 0.2, 1.0);
                gl.clear(gl.COLOR_BUFFER_BIT | gl.DEPTH_BUFFER_BIT);

                // Enable depth testing
                gl.enable(gl.DEPTH_TEST);
                gl.depthFunc(gl.LEQUAL);

                // Draw the cube using indices
                gl.drawElements(gl.TRIANGLES, data.indices.length, gl.UNSIGNED_SHORT, 0);

                statusMessage.textContent = 'Rendering complete!';

            } catch (error) {
                console.error('Error fetching data or rendering:', error);
                statusMessage.textContent = `Error: ${error.message}`;
            }

            // A simple mat4 library for matrix operations (needed for WebGL)
            const mat4 = {
                create: function() { return new Float32Array(16); },
                identity: function(out) {
                    out[0] = 1; out[1] = 0; out[2] = 0; out[3] = 0;
                    out[4] = 0; out[5] = 1; out[6] = 0; out[7] = 0;
                    out[8] = 0; out[9] = 0; out[10] = 1; out[11] = 0;
                    out[12] = 0; out[13] = 0; out[14] = 0; out[15] = 1;
                    return out;
                },
                translate: function(out, a, v) {
                    const x = v[0], y = v[1], z = v[2];
                    out[0] = a[0]; out[1] = a[1]; out[2] = a[2]; out[3] = a[3];
                    out[4] = a[4]; out[5] = a[5]; out[6] = a[6]; out[7] = a[7];
                    out[8] = a[8]; out[9] = a[9]; out[10] = a[10]; out[11] = a[11];
                    let a0 = a[0], a1 = a[1], a2 = a[2], a3 = a[3];
                    let a4 = a[4], a5 = a[5], a6 = a[6], a7 = a[7];
                    let a8 = a[8], a9 = a[9], a10 = a[10], a11 = a[11];
                    out[12] = a0 * x + a4 * y + a8 * z + a[12];
                    out[13] = a1 * x + a5 * y + a9 * z + a[13];
                    out[14] = a2 * x + a6 * y + a10 * z + a[14];
                    out[15] = a3 * x + a7 * y + a11 * z + a[15];
                    return out;
                },
                perspective: function(out, fov, aspect, near, far) {
                    const f = 1.0 / Math.tan(fov / 2);
                    out[0] = f / aspect;
                    out[1] = 0;
                    out[2] = 0;
                    out[3] = 0;
                    out[4] = 0;
                    out[5] = f;
                    out[6] = 0;
                    out[7] = 0;
                    out[8] = 0;
                    out[9] = 0;
                    out[10] = (near + far) / (near - far);
                    out[11] = -1;
                    out[12] = 0;
                    out[13] = 0;
                    out[14] = (2 * near * far) / (near - far);
                    out[15] = 0;
                    return out;
                }
            };

        });
</script>