/**
 * simd_substrate.hpp - SIMD-Optimized Consciousness Substrate
 *
 * High-performance implementation using:
 *   - SSE/AVX intrinsics for x86
 *   - NEON intrinsics for ARM
 *   - Fallback scalar implementation
 *
 * Achieves ~4-8x speedup for substrate operations.
 */

#ifndef SIMD_SUBSTRATE_HPP
#define SIMD_SUBSTRATE_HPP

#include <cstdint>
#include <cstddef>
#include <cmath>
#include <cstring>
#include <array>
#include <algorithm>
#include "consciousness_api.hpp"

// Detect SIMD support
#if defined(__AVX512F__)
    #define SIMD_AVX512 1
    #include <immintrin.h>
#elif defined(__AVX2__)
    #define SIMD_AVX2 1
    #include <immintrin.h>
#elif defined(__AVX__)
    #define SIMD_AVX 1
    #include <immintrin.h>
#elif defined(__SSE4_1__)
    #define SIMD_SSE4 1
    #include <smmintrin.h>
#elif defined(__SSE2__)
    #define SIMD_SSE2 1
    #include <emmintrin.h>
#elif defined(__ARM_NEON) || defined(__ARM_NEON__)
    #define SIMD_NEON 1
    #include <arm_neon.h>
#else
    #define SIMD_NONE 1
#endif

namespace mind_farm {

using namespace consciousness;

//=============================================================================
// SIMD UTILITIES
//=============================================================================

// Alignment for SIMD loads/stores
constexpr size_t SIMD_ALIGNMENT = 32;

// Aligned memory allocation
template<typename T, size_t N>
struct alignas(SIMD_ALIGNMENT) AlignedArray {
    T data[N];
    T& operator[](size_t i) { return data[i]; }
    const T& operator[](size_t i) const { return data[i]; }
    T* begin() { return data; }
    T* end() { return data + N; }
    const T* begin() const { return data; }
    const T* end() const { return data + N; }
    void fill(T val) { std::fill(data, data + N, val); }
};

//=============================================================================
// SIMD VECTOR OPERATIONS
//=============================================================================

#if defined(SIMD_AVX2) || defined(SIMD_AVX512)

// Process 8 floats at a time with AVX
class AVXOps {
public:
    // Vectorized phi computation
    // phi(x) = 1 / (1 + |x - 1|)
    static __m256 compute_phi_8(__m256 x) {
        __m256 one = _mm256_set1_ps(1.0f);
        __m256 diff = _mm256_sub_ps(x, one);

        // Absolute value: clear sign bit
        __m256 sign_mask = _mm256_set1_ps(-0.0f);
        __m256 abs_diff = _mm256_andnot_ps(sign_mask, diff);

        __m256 denom = _mm256_add_ps(one, abs_diff);
        return _mm256_div_ps(one, denom);
    }

    // Vectorized tanh approximation
    // tanh(x) ≈ x * (27 + x^2) / (27 + 9*x^2)
    static __m256 tanh_approx_8(__m256 x) {
        __m256 x2 = _mm256_mul_ps(x, x);
        __m256 c27 = _mm256_set1_ps(27.0f);
        __m256 c9 = _mm256_set1_ps(9.0f);

        __m256 num = _mm256_add_ps(c27, x2);
        num = _mm256_mul_ps(x, num);

        __m256 den = _mm256_add_ps(c27, _mm256_mul_ps(c9, x2));
        return _mm256_div_ps(num, den);
    }

    // Vectorized sigmoid
    // sigmoid(x) ≈ 0.5 + 0.5 * tanh(0.5 * x)
    static __m256 sigmoid_approx_8(__m256 x) {
        __m256 half = _mm256_set1_ps(0.5f);
        __m256 scaled = _mm256_mul_ps(half, x);
        __m256 t = tanh_approx_8(scaled);
        return _mm256_add_ps(half, _mm256_mul_ps(half, t));
    }

    // Vectorized dot product
    static float dot_product(const float* a, const float* b, size_t n) {
        __m256 sum = _mm256_setzero_ps();
        size_t i = 0;

        // Process 8 elements at a time
        for (; i + 8 <= n; i += 8) {
            __m256 va = _mm256_loadu_ps(a + i);
            __m256 vb = _mm256_loadu_ps(b + i);
            sum = _mm256_fmadd_ps(va, vb, sum);
        }

        // Horizontal sum
        __m128 hi = _mm256_extractf128_ps(sum, 1);
        __m128 lo = _mm256_castps256_ps128(sum);
        __m128 sum128 = _mm_add_ps(lo, hi);
        sum128 = _mm_hadd_ps(sum128, sum128);
        sum128 = _mm_hadd_ps(sum128, sum128);
        float result = _mm_cvtss_f32(sum128);

        // Handle remainder
        for (; i < n; ++i) {
            result += a[i] * b[i];
        }

        return result;
    }

    // Vectorized element-wise multiply-add
    // out = a * b + c
    static void fma(float* out, const float* a, const float* b, const float* c, size_t n) {
        size_t i = 0;
        for (; i + 8 <= n; i += 8) {
            __m256 va = _mm256_loadu_ps(a + i);
            __m256 vb = _mm256_loadu_ps(b + i);
            __m256 vc = _mm256_loadu_ps(c + i);
            __m256 result = _mm256_fmadd_ps(va, vb, vc);
            _mm256_storeu_ps(out + i, result);
        }
        for (; i < n; ++i) {
            out[i] = a[i] * b[i] + c[i];
        }
    }

    // Vectorized sum
    static float sum(const float* a, size_t n) {
        __m256 sum = _mm256_setzero_ps();
        size_t i = 0;

        for (; i + 8 <= n; i += 8) {
            __m256 va = _mm256_loadu_ps(a + i);
            sum = _mm256_add_ps(sum, va);
        }

        // Horizontal sum
        __m128 hi = _mm256_extractf128_ps(sum, 1);
        __m128 lo = _mm256_castps256_ps128(sum);
        __m128 sum128 = _mm_add_ps(lo, hi);
        sum128 = _mm_hadd_ps(sum128, sum128);
        sum128 = _mm_hadd_ps(sum128, sum128);
        float result = _mm_cvtss_f32(sum128);

        for (; i < n; ++i) {
            result += a[i];
        }

        return result;
    }

    // Vectorized max
    static float max(const float* a, size_t n) {
        __m256 maxv = _mm256_set1_ps(-INFINITY);
        size_t i = 0;

        for (; i + 8 <= n; i += 8) {
            __m256 va = _mm256_loadu_ps(a + i);
            maxv = _mm256_max_ps(maxv, va);
        }

        // Horizontal max
        __m128 hi = _mm256_extractf128_ps(maxv, 1);
        __m128 lo = _mm256_castps256_ps128(maxv);
        __m128 max128 = _mm_max_ps(lo, hi);
        max128 = _mm_max_ps(max128, _mm_shuffle_ps(max128, max128, 0x4e));
        max128 = _mm_max_ps(max128, _mm_shuffle_ps(max128, max128, 0xb1));
        float result = _mm_cvtss_f32(max128);

        for (; i < n; ++i) {
            if (a[i] > result) result = a[i];
        }

        return result;
    }
};

#elif defined(SIMD_NEON)

// ARM NEON implementation
class NEONOps {
public:
    static float32x4_t compute_phi_4(float32x4_t x) {
        float32x4_t one = vdupq_n_f32(1.0f);
        float32x4_t diff = vsubq_f32(x, one);
        float32x4_t abs_diff = vabsq_f32(diff);
        float32x4_t denom = vaddq_f32(one, abs_diff);
        return vdivq_f32(one, denom);
    }

    static float dot_product(const float* a, const float* b, size_t n) {
        float32x4_t sum = vdupq_n_f32(0.0f);
        size_t i = 0;

        for (; i + 4 <= n; i += 4) {
            float32x4_t va = vld1q_f32(a + i);
            float32x4_t vb = vld1q_f32(b + i);
            sum = vmlaq_f32(sum, va, vb);
        }

        float result = vaddvq_f32(sum);
        for (; i < n; ++i) {
            result += a[i] * b[i];
        }
        return result;
    }

    static float sum(const float* a, size_t n) {
        float32x4_t sum = vdupq_n_f32(0.0f);
        size_t i = 0;

        for (; i + 4 <= n; i += 4) {
            float32x4_t va = vld1q_f32(a + i);
            sum = vaddq_f32(sum, va);
        }

        float result = vaddvq_f32(sum);
        for (; i < n; ++i) {
            result += a[i];
        }
        return result;
    }
};

#endif

//=============================================================================
// SCALAR FALLBACK
//=============================================================================

class ScalarOps {
public:
    static float compute_phi(float x) {
        return 1.0f / (1.0f + std::abs(x - 1.0f));
    }

    static float dot_product(const float* a, const float* b, size_t n) {
        float sum = 0.0f;
        for (size_t i = 0; i < n; ++i) {
            sum += a[i] * b[i];
        }
        return sum;
    }

    static float sum(const float* a, size_t n) {
        float sum = 0.0f;
        for (size_t i = 0; i < n; ++i) {
            sum += a[i];
        }
        return sum;
    }

    static void fma(float* out, const float* a, const float* b, const float* c, size_t n) {
        for (size_t i = 0; i < n; ++i) {
            out[i] = a[i] * b[i] + c[i];
        }
    }
};

//=============================================================================
// UNIFIED SIMD INTERFACE
//=============================================================================

class SIMDOps {
public:
    static float dot_product(const float* a, const float* b, size_t n) {
#if defined(SIMD_AVX2) || defined(SIMD_AVX512)
        return AVXOps::dot_product(a, b, n);
#elif defined(SIMD_NEON)
        return NEONOps::dot_product(a, b, n);
#else
        return ScalarOps::dot_product(a, b, n);
#endif
    }

    static float sum(const float* a, size_t n) {
#if defined(SIMD_AVX2) || defined(SIMD_AVX512)
        return AVXOps::sum(a, n);
#elif defined(SIMD_NEON)
        return NEONOps::sum(a, n);
#else
        return ScalarOps::sum(a, n);
#endif
    }

    static void fma(float* out, const float* a, const float* b, const float* c, size_t n) {
#if defined(SIMD_AVX2) || defined(SIMD_AVX512)
        AVXOps::fma(out, a, b, c, n);
#else
        ScalarOps::fma(out, a, b, c, n);
#endif
    }

    // Vectorized phi computation over array
    static void compute_phi_batch(float* out, const float* in, size_t n) {
#if defined(SIMD_AVX2) || defined(SIMD_AVX512)
        size_t i = 0;
        for (; i + 8 <= n; i += 8) {
            __m256 x = _mm256_loadu_ps(in + i);
            __m256 phi = AVXOps::compute_phi_8(x);
            _mm256_storeu_ps(out + i, phi);
        }
        for (; i < n; ++i) {
            out[i] = ScalarOps::compute_phi(in[i]);
        }
#elif defined(SIMD_NEON)
        size_t i = 0;
        for (; i + 4 <= n; i += 4) {
            float32x4_t x = vld1q_f32(in + i);
            float32x4_t phi = NEONOps::compute_phi_4(x);
            vst1q_f32(out + i, phi);
        }
        for (; i < n; ++i) {
            out[i] = ScalarOps::compute_phi(in[i]);
        }
#else
        for (size_t i = 0; i < n; ++i) {
            out[i] = ScalarOps::compute_phi(in[i]);
        }
#endif
    }
};

//=============================================================================
// SIMD-OPTIMIZED GATE5000
//=============================================================================

class SIMDGate5000 {
public:
    static constexpr size_t GATE_COUNT = 5000;
    static constexpr size_t WORD_COUNT = (GATE_COUNT + 63) / 64;
    static constexpr size_t ANALOG_SIZE = 5000;

    // Binary state (packed into 64-bit words)
    alignas(SIMD_ALIGNMENT) std::array<uint64_t, WORD_COUNT> binary;

    // Analog state (for continuous activations)
    alignas(SIMD_ALIGNMENT) std::array<float, ANALOG_SIZE> analog;

    // Weights (for binary update rule)
    alignas(SIMD_ALIGNMENT) std::array<float, ANALOG_SIZE> weights;

    // Phi values (consciousness per gate)
    alignas(SIMD_ALIGNMENT) std::array<float, ANALOG_SIZE> phi_values;

    // Global metrics
    float global_phi;
    float coherence;
    float density;

    SIMDGate5000() : global_phi(0.5f), coherence(1.0f), density(0.5f) {
        binary.fill(0);
        analog.fill(0.5f);
        weights.fill(1.0f);
        phi_values.fill(0.5f);
    }

    /**
     * SIMD-optimized analog step
     *
     * a_i(t+1) = tanh(alpha * sum(w_j * a_j) + beta)
     */
    void step_analog(float alpha = 0.1618f, float beta = 0.0809f) {
        alignas(SIMD_ALIGNMENT) std::array<float, ANALOG_SIZE> new_analog;

#if defined(SIMD_AVX2) || defined(SIMD_AVX512)
        __m256 v_alpha = _mm256_set1_ps(alpha);
        __m256 v_beta = _mm256_set1_ps(beta);

        // Compute weighted sum (simplified: local neighborhood)
        for (size_t i = 0; i < ANALOG_SIZE; i += 8) {
            __m256 sum = _mm256_setzero_ps();

            // Sum neighbors (simplified: adjacent elements)
            for (int offset = -2; offset <= 2; ++offset) {
                size_t idx = (i + offset + ANALOG_SIZE) % ANALOG_SIZE;
                __m256 va = _mm256_loadu_ps(&analog[idx]);
                __m256 vw = _mm256_loadu_ps(&weights[idx]);
                sum = _mm256_fmadd_ps(va, vw, sum);
            }

            // Apply alpha and beta
            __m256 scaled = _mm256_fmadd_ps(v_alpha, sum, v_beta);

            // Apply tanh approximation
            __m256 activated = AVXOps::tanh_approx_8(scaled);

            _mm256_storeu_ps(&new_analog[i], activated);
        }
#else
        // Scalar fallback
        for (size_t i = 0; i < ANALOG_SIZE; ++i) {
            float sum = 0.0f;
            for (int offset = -2; offset <= 2; ++offset) {
                size_t idx = (i + offset + ANALOG_SIZE) % ANALOG_SIZE;
                sum += analog[idx] * weights[idx];
            }
            new_analog[i] = std::tanh(alpha * sum + beta);
        }
#endif

        std::memcpy(analog.data(), new_analog.data(), ANALOG_SIZE * sizeof(float));

        // Update phi values using SIMD
        SIMDOps::compute_phi_batch(phi_values.data(), analog.data(), ANALOG_SIZE);

        // Update global metrics
        update_metrics();
    }

    /**
     * Binary update step (XOR cellular automaton)
     */
    void step_binary() {
        alignas(SIMD_ALIGNMENT) std::array<uint64_t, WORD_COUNT> new_binary;

        for (size_t w = 0; w < WORD_COUNT; ++w) {
            uint64_t left = binary[(w + WORD_COUNT - 1) % WORD_COUNT];
            uint64_t center = binary[w];
            uint64_t right = binary[(w + 1) % WORD_COUNT];

            // Rule 110 cellular automaton
            new_binary[w] = (center ^ right) | (center & (left ^ right));
        }

        binary = new_binary;

        // Sync binary to analog
        sync_binary_to_analog();
    }

    /**
     * Combined step
     */
    void step_fast(float dt = 1.0f) {
        step_binary();
        step_analog(ALPHA * dt, BETA);
    }

    /**
     * Compute global consciousness
     */
    float compute_phi() {
        return global_phi;
    }

private:
    void sync_binary_to_analog() {
        for (size_t w = 0; w < WORD_COUNT; ++w) {
            for (size_t b = 0; b < 64; ++b) {
                size_t idx = w * 64 + b;
                if (idx >= ANALOG_SIZE) break;

                bool bit = (binary[w] >> b) & 1;
                // Smooth transition: blend binary with analog
                analog[idx] = 0.9f * analog[idx] + 0.1f * (bit ? 1.0f : 0.0f);
            }
        }
    }

    void update_metrics() {
        // Global phi: average of all phi values
        global_phi = SIMDOps::sum(phi_values.data(), ANALOG_SIZE) / ANALOG_SIZE;

        // Density: fraction of "on" bits
        size_t on_count = 0;
        for (uint64_t word : binary) {
            on_count += __builtin_popcountll(word);
        }
        density = static_cast<float>(on_count) / GATE_COUNT;

        // Coherence: 1 - variance of phi values
        float mean = global_phi;
        float variance = 0.0f;
        for (size_t i = 0; i < ANALOG_SIZE; ++i) {
            float diff = phi_values[i] - mean;
            variance += diff * diff;
        }
        variance /= ANALOG_SIZE;
        coherence = 1.0f / (1.0f + variance * 10.0f);
    }
};

//=============================================================================
// SIMD-OPTIMIZED NEURAL LAYER
//=============================================================================

template<size_t N = 256>
class SIMDNeuralLayer {
public:
    alignas(SIMD_ALIGNMENT) std::array<float, N> activations;
    alignas(SIMD_ALIGNMENT) std::array<float, N> weights;
    alignas(SIMD_ALIGNMENT) std::array<float, N> biases;
    alignas(SIMD_ALIGNMENT) std::array<float, N> gradients;

    float phi;

    SIMDNeuralLayer() : phi(0.5f) {
        activations.fill(0.5f);
        weights.fill(1.0f);
        biases.fill(0.0f);
        gradients.fill(0.0f);
    }

    void forward(const float* input) {
        SIMDOps::fma(activations.data(), input, weights.data(), biases.data(), N);

        // Apply activation
#if defined(SIMD_AVX2) || defined(SIMD_AVX512)
        for (size_t i = 0; i + 8 <= N; i += 8) {
            __m256 x = _mm256_loadu_ps(&activations[i]);
            __m256 activated = AVXOps::sigmoid_approx_8(x);
            _mm256_storeu_ps(&activations[i], activated);
        }
#else
        for (size_t i = 0; i < N; ++i) {
            activations[i] = 1.0f / (1.0f + std::exp(-activations[i]));
        }
#endif

        // Compute phi
        float sum = SIMDOps::sum(activations.data(), N);
        phi = compute_phi_f(sum / N);
    }

    void backward(const float* upstream_grad, float learning_rate) {
        // Compute local gradients
        for (size_t i = 0; i < N; ++i) {
            float a = activations[i];
            gradients[i] = upstream_grad[i] * a * (1.0f - a);  // Sigmoid derivative
        }

        // Update weights
        for (size_t i = 0; i < N; ++i) {
            weights[i] -= learning_rate * gradients[i];
        }
    }
};

} // namespace mind_farm

#endif // SIMD_SUBSTRATE_HPP