rgdl_engine.py

"""
Universal Binary Principle (UBP) Framework v3.1.1 - Enhanced RGDL Engine Module
Author: Euan Craig, New Zealand
Date: 18 August 2025

This module implements the Resonance Geometry Definition Language (RGDL)
Geometric Execution Engine, providing dynamic geometry generation through
emergent behavior of binary toggles operating under specific resonance
frequencies and coherence constraints.

Enhanced for v3.1.1 with improved integration with v3.1 components and
better performance optimization.
"""

import numpy as np
from typing import Dict, Any, List, Tuple, Optional, Union, Callable
from dataclasses import dataclass
import json
import math
import time
from scipy.spatial import ConvexHull, Voronoi
from scipy.spatial.distance import pdist, squareform
from scipy.optimize import minimize

try:
    from .system_constants import UBPConstants # Corrected import from 'core' to 'system_constants'
    from .bits import Bitfield, OffBit # Corrected import from 'bitfield' to 'bits'
    from .toggle_algebra import ToggleAlgebra
    from .hex_dictionary import HexDictionary
except ImportError:
    from system_constants import UBPConstants # Corrected import from 'core' to 'system_constants'
    from bits import Bitfield, OffBit # Corrected import from 'bitfield' to 'bits'
    from toggle_algebra import ToggleAlgebra
    from hex_dictionary import HexDictionary


@dataclass
class GeometricPrimitive:
    """A geometric primitive generated by RGDL."""
    primitive_type: str
    coordinates: np.ndarray
    properties: Dict[str, Any]
    resonance_frequency: float
    coherence_level: float
    generation_method: str
    stability_score: float
    creation_timestamp: float
    nrci_score: float = 0.0


@dataclass
class RGDLMetrics:
    """Performance and quality metrics for RGDL operations."""
    total_primitives_generated: int
    average_coherence: float
    average_stability: float
    geometric_complexity: float
    resonance_distribution: Dict[str, int]
    generation_time: float
    memory_usage_mb: float
    nrci_average: float = 0.0


@dataclass
class GeometricField:
    """A field of geometric primitives with spatial relationships."""
    field_name: str
    primitives: List[GeometricPrimitive]
    spatial_bounds: Tuple[Tuple[float, float], Tuple[float, float], Tuple[float, float]]
    field_coherence: float
    resonance_pattern: np.ndarray
    interaction_matrix: np.ndarray
    field_energy: float = 0.0


class RGDLEngine:
    """
    Enhanced Resonance Geometry Definition Language (RGDL) Execution Engine.
    
    This engine generates geometric primitives through the emergent behavior
    of binary toggles operating under specific resonance frequencies and
    coherence constraints within the UBP framework.
    
    Enhanced for v3.1 with better integration and performance.
    """
    
    def __init__(self, bitfield_instance: Optional[Bitfield] = None,
                 toggle_algebra_instance: Optional[ToggleAlgebra] = None,
                 hex_dictionary_instance: Optional[HexDictionary] = None):
        """
        Initialize the RGDL Engine.
        
        Args:
            bitfield_instance: Optional Bitfield instance for geometric operations
            toggle_algebra_instance: Optional ToggleAlgebra instance for computations
            hex_dictionary_instance: Optional HexDictionary for data storage
        """
        self.bitfield = bitfield_instance
        self.toggle_algebra = toggle_algebra_instance
        self.hex_dictionary = hex_dictionary_instance or HexDictionary()
        
        # Geometric primitive generators
        self.primitive_generators = {
            'point': self._generate_point,
            'line': self._generate_line,
            'triangle': self._generate_triangle,
            'tetrahedron': self._generate_tetrahedron,
            'cube': self._generate_cube,
            'sphere': self._generate_sphere,
            'torus': self._generate_torus,
            'fractal': self._generate_fractal,
            'resonance_surface': self._generate_resonance_surface,
            'quantum_geometry': self._generate_quantum_geometry,  # New for v3.1
            'htr_structure': self._generate_htr_structure,  # New for v3.1
            'crv_manifold': self._generate_crv_manifold,  # New for v3.1
        }
        
        # Resonance frequency mappings
        # Assuming UBPConstants provides these values, or using hardcoded examples for standalone testing
        self.resonance_frequencies = {
            'quantum': getattr(UBPConstants, 'CRV_QUANTUM_BASE', 0.2265234857),
            'electromagnetic': getattr(UBPConstants, 'CRV_ELECTROMAGNETIC_BASE', 3.141593),
            'gravitational': getattr(UBPConstants, 'CRV_GRAVITATIONAL_BASE', 100.0),
            'biological': getattr(UBPConstants, 'CRV_BIOLOGICAL_BASE', 10.0),
            'cosmological': getattr(UBPConstants, 'CRV_COSMOLOGICAL_BASE', 0.83203682),
            'nuclear': getattr(UBPConstants, 'CRV_NUCLEAR_BASE', 1.2356e20),
            'optical': getattr(UBPConstants, 'CRV_OPTICAL_BASE', 5.0e14)
        }
        
        # Generated primitives storage
        self.primitives: List[GeometricPrimitive] = []
        self.geometric_fields: Dict[str, GeometricField] = {}
        
        # Performance metrics
        self.metrics = RGDLMetrics(
            total_primitives_generated=0,
            average_coherence=0.0,
            average_stability=0.0,
            geometric_complexity=0.0,
            resonance_distribution={},
            generation_time=0.0,
            memory_usage_mb=0.0
        )
        
        # Geometry cache for performance
        self.geometry_cache: Dict[str, GeometricPrimitive] = {}
        
        print("✅ RGDL Engine v3.1 Initialized")
        print(f"   Supported Primitives: {len(self.primitive_generators)}")
        print(f"   Resonance Frequencies: {len(self.resonance_frequencies)}")
        print(f"   HexDictionary Integration: {'Enabled' if self.hex_dictionary else 'Disabled'}")
    
    # ========================================================================
    # CORE GEOMETRY GENERATION METHODS
    # ========================================================================
    
    def generate_primitive(self, primitive_type: str, 
                          resonance_realm: str = 'electromagnetic',
                          parameters: Optional[Dict[str, Any]] = None) -> GeometricPrimitive:
        """
        Generate a geometric primitive of the specified type.
        
        Args:
            primitive_type: Type of primitive to generate
            resonance_realm: Realm for resonance frequency selection
            parameters: Optional parameters for generation
            
        Returns:
            Generated GeometricPrimitive
        """
        start_time = time.time()
        
        if primitive_type not in self.primitive_generators:
            raise ValueError(f"Unsupported primitive type: {primitive_type}")
        
        if resonance_realm not in self.resonance_frequencies:
            raise ValueError(f"Unsupported resonance realm: {resonance_realm}")
        
        # Get resonance frequency
        resonance_freq = self.resonance_frequencies[resonance_realm]
        
        # Generate cache key
        cache_key = f"{primitive_type}_{resonance_realm}_{hash(str(parameters))}"
        
        # Check cache first
        if cache_key in self.geometry_cache:
            cached_primitive = self.geometry_cache[cache_key]
            cached_primitive.creation_timestamp = time.time()
            return cached_primitive
        
        # Generate the primitive
        generator = self.primitive_generators[primitive_type]
        primitive = generator(resonance_freq, parameters or {})
        
        # Calculate coherence and stability
        primitive.coherence_level = self._calculate_coherence(primitive)
        primitive.stability_score = self._calculate_stability(primitive)
        primitive.nrci_score = self._calculate_nrci(primitive)
        primitive.creation_timestamp = time.time()
        
        # Store in cache
        self.geometry_cache[cache_key] = primitive
        
        # Store in HexDictionary if available
        if self.hex_dictionary:
            geometry_data = {
                'type': primitive_type,
                'coordinates': primitive.coordinates.tolist(),
                'properties': primitive.properties,
                'resonance_frequency': primitive.resonance_frequency,
                'coherence_level': primitive.coherence_level,
                'stability_score': primitive.stability_score,
                'nrci_score': primitive.nrci_score
            }
            self.hex_dictionary.store(geometry_data, 'json', 
                                    {'primitive_type': primitive_type, 
                                     'resonance_realm': resonance_realm})
        
        # Update metrics
        self.primitives.append(primitive)
        self._update_metrics(primitive, time.time() - start_time)
        
        return primitive
    
    def _generate_point(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a resonance-influenced point."""
        # Use resonance frequency to influence position
        phase = params.get('phase', 0.0)
        amplitude = params.get('amplitude', 1.0)
        
        x = amplitude * np.cos(2 * np.pi * resonance_freq * phase)
        y = amplitude * np.sin(2 * np.pi * resonance_freq * phase)
        z = amplitude * np.cos(np.pi * resonance_freq * phase)
        
        coordinates = np.array([x, y, z])
        
        return GeometricPrimitive(
            primitive_type='point',
            coordinates=coordinates,
            properties={'amplitude': amplitude, 'phase': phase},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,  # Will be calculated later
            generation_method='resonance_oscillation',
            stability_score=0.0,  # Will be calculated later
            creation_timestamp=time.time()
        )
    
    def _generate_line(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a resonance-influenced line segment."""
        length = params.get('length', 1.0)
        direction = params.get('direction', np.array([1, 0, 0]))
        start_point = params.get('start_point', np.array([0, 0, 0]))
        
        # Normalize direction
        direction = direction / np.linalg.norm(direction)
        
        # Create line points influenced by resonance
        num_points = params.get('num_points', 100)
        t_values = np.linspace(0, length, num_points)
        
        # Add resonance-based perturbation
        perturbation_amplitude = params.get('perturbation', 0.01)
        perturbation = perturbation_amplitude * np.sin(2 * np.pi * resonance_freq * t_values)
        
        coordinates = np.array([start_point + t * direction + 
                              perturbation[i] * np.array([0, 1, 0]) 
                              for i, t in enumerate(t_values)])
        
        return GeometricPrimitive(
            primitive_type='line',
            coordinates=coordinates,
            properties={'length': length, 'direction': direction.tolist(), 
                       'perturbation_amplitude': perturbation_amplitude},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='resonance_perturbation',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_triangle(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a resonance-influenced triangle."""
        center = params.get('center', np.array([0, 0, 0]))
        radius = params.get('radius', 1.0)
        
        # Generate triangle vertices with resonance influence
        angles = np.array([0, 2*np.pi/3, 4*np.pi/3])
        resonance_modulation = 1 + 0.1 * np.sin(2 * np.pi * resonance_freq * angles)
        
        vertices = []
        for i, angle in enumerate(angles):
            x = center[0] + radius * resonance_modulation[i] * np.cos(angle)
            y = center[1] + radius * resonance_modulation[i] * np.sin(angle)
            z = center[2]
            vertices.append([x, y, z])
        
        coordinates = np.array(vertices)
        
        return GeometricPrimitive(
            primitive_type='triangle',
            coordinates=coordinates,
            properties={'center': center.tolist(), 'radius': radius, 
                       'resonance_modulation': resonance_modulation.tolist()},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='resonance_modulation',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_tetrahedron(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a resonance-influenced tetrahedron."""
        center = params.get('center', np.array([0, 0, 0]))
        edge_length = params.get('edge_length', 1.0)
        
        # Standard tetrahedron vertices
        a = edge_length / np.sqrt(2)
        vertices = np.array([
            [a, a, a],
            [a, -a, -a],
            [-a, a, -a],
            [-a, -a, a]
        ]) + center
        
        # Apply resonance-based deformation
        deformation_factor = 0.1 * np.sin(2 * np.pi * resonance_freq)
        vertices *= (1 + deformation_factor)
        
        return GeometricPrimitive(
            primitive_type='tetrahedron',
            coordinates=vertices,
            properties={'center': center.tolist(), 'edge_length': edge_length,
                       'deformation_factor': deformation_factor},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='resonance_deformation',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_cube(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a resonance-influenced cube."""
        center = params.get('center', np.array([0, 0, 0]))
        side_length = params.get('side_length', 1.0)
        
        # Generate cube vertices
        half_side = side_length / 2
        vertices = []
        
        for x in [-half_side, half_side]:
            for y in [-half_side, half_side]:
                for z in [-half_side, half_side]:
                    # Apply resonance-based position adjustment
                    resonance_shift = 0.05 * np.sin(2 * np.pi * resonance_freq * (x + y + z))
                    vertex = center + np.array([x, y, z]) * (1 + resonance_shift)
                    vertices.append(vertex)
        
        coordinates = np.array(vertices)
        
        return GeometricPrimitive(
            primitive_type='cube',
            coordinates=coordinates,
            properties={'center': center.tolist(), 'side_length': side_length},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='resonance_vertex_shift',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_sphere(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a resonance-influenced sphere."""
        center = params.get('center', np.array([0, 0, 0]))
        radius = params.get('radius', 1.0)
        resolution = params.get('resolution', 50)
        
        # Generate sphere points
        phi = np.linspace(0, np.pi, resolution)
        theta = np.linspace(0, 2*np.pi, resolution)
        
        vertices = []
        for p in phi:
            for t in theta:
                # Apply resonance-based radius modulation
                r_modulated = radius * (1 + 0.1 * np.sin(2 * np.pi * resonance_freq * (p + t)))
                
                x = center[0] + r_modulated * np.sin(p) * np.cos(t)
                y = center[1] + r_modulated * np.sin(p) * np.sin(t)
                z = center[2] + r_modulated * np.cos(p)
                
                vertices.append([x, y, z])
        
        coordinates = np.array(vertices)
        
        return GeometricPrimitive(
            primitive_type='sphere',
            coordinates=coordinates,
            properties={'center': center.tolist(), 'radius': radius, 'resolution': resolution},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='resonance_radius_modulation',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_torus(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a resonance-influenced torus."""
        center = params.get('center', np.array([0, 0, 0]))
        major_radius = params.get('major_radius', 1.0)
        minor_radius = params.get('minor_radius', 0.3)
        resolution = params.get('resolution', 30)
        
        # Generate torus points
        u = np.linspace(0, 2*np.pi, resolution)
        v = np.linspace(0, 2*np.pi, resolution)
        
        vertices = []
        for u_val in u:
            for v_val in v:
                # Apply resonance-based modulation
                resonance_factor = 1 + 0.1 * np.sin(2 * np.pi * resonance_freq * (u_val + v_val))
                
                x = center[0] + (major_radius + minor_radius * np.cos(v_val)) * np.cos(u_val) * resonance_factor
                y = center[1] + (major_radius + minor_radius * np.cos(v_val)) * np.sin(u_val) * resonance_factor
                z = center[2] + minor_radius * np.sin(v_val) * resonance_factor
                
                vertices.append([x, y, z])
        
        coordinates = np.array(vertices)
        
        return GeometricPrimitive(
            primitive_type='torus',
            coordinates=coordinates,
            properties={'center': center.tolist(), 'major_radius': major_radius, 
                       'minor_radius': minor_radius, 'resolution': resolution},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='resonance_torus_modulation',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_fractal(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a resonance-influenced fractal structure."""
        iterations = params.get('iterations', 5)
        scale_factor = params.get('scale_factor', 0.5)
        base_shape = params.get('base_shape', 'triangle')
        
        # Start with base shape
        if base_shape == 'triangle':
            base_vertices = np.array([[0, 0, 0], [1, 0, 0], [0.5, np.sqrt(3)/2, 0]])
        else:
            base_vertices = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]])
        
        # Apply fractal generation with resonance influence
        all_vertices = [base_vertices]
        
        for iteration in range(iterations):
            new_vertices = []
            resonance_scale = 1 + 0.1 * np.sin(2 * np.pi * resonance_freq * iteration)
            
            for vertex_set in all_vertices:
                # Create smaller copies at each vertex
                for vertex in vertex_set:
                    scaled_shape = vertex + (vertex_set - vertex_set[0]) * scale_factor * resonance_scale
                    new_vertices.append(scaled_shape)
            
            all_vertices.extend(new_vertices)
        
        # Flatten all vertices
        coordinates = np.vstack(all_vertices)
        
        return GeometricPrimitive(
            primitive_type='fractal',
            coordinates=coordinates,
            properties={'iterations': iterations, 'scale_factor': scale_factor, 
                       'base_shape': base_shape},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='resonance_fractal_scaling',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_resonance_surface(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate a surface based on resonance patterns."""
        x_range = params.get('x_range', (-2, 2))
        y_range = params.get('y_range', (-2, 2))
        resolution = params.get('resolution', 50)
        
        x = np.linspace(x_range[0], x_range[1], resolution)
        y = np.linspace(y_range[0], y_range[1], resolution)
        X, Y = np.meshgrid(x, y)
        
        # Generate resonance-based surface
        Z = np.sin(2 * np.pi * resonance_freq * X) * np.cos(2 * np.pi * resonance_freq * Y)
        Z += 0.5 * np.sin(4 * np.pi * resonance_freq * np.sqrt(X**2 + Y**2))
        
        # Convert to vertex array
        vertices = []
        for i in range(resolution):
            for j in range(resolution):
                vertices.append([X[i, j], Y[i, j], Z[i, j]])
        
        coordinates = np.array(vertices)
        
        return GeometricPrimitive(
            primitive_type='resonance_surface',
            coordinates=coordinates,
            properties={'x_range': x_range, 'y_range': y_range, 'resolution': resolution},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='resonance_wave_interference',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    # New v3.1 generators
    def _generate_quantum_geometry(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate quantum-influenced geometry using UBP principles."""
        quantum_states = params.get('quantum_states', 4)
        superposition_factor = params.get('superposition_factor', 0.5)
        
        # Generate quantum state positions
        vertices = []
        for state in range(quantum_states):
            # Use quantum CRV for positioning
            phase = 2 * np.pi * state / quantum_states
            # Assuming UBPConstants.CRV_QUANTUM_BASE exists or fallback
            quantum_crv = getattr(UBPConstants, 'CRV_QUANTUM_BASE', 0.2265234857)
            
            x = np.cos(phase) * (1 + superposition_factor * np.sin(2 * np.pi * quantum_crv * state))
            y = np.sin(phase) * (1 + superposition_factor * np.cos(2 * np.pi * quantum_crv * state))
            z = superposition_factor * np.sin(2 * np.pi * quantum_crv * phase)
            
            vertices.append([x, y, z])
        
        coordinates = np.array(vertices)
        
        return GeometricPrimitive(
            primitive_type='quantum_geometry',
            coordinates=coordinates,
            properties={'quantum_states': quantum_states, 'superposition_factor': superposition_factor},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='quantum_superposition',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_htr_structure(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """
        Generate HTR (Harmonic Toggle Resonance) influenced structure.
        The geometry is now explicitly influenced by the NRCI value from HTR Engine.
        """
        harmonic_order = params.get('harmonic_order', 3)
        nrci_value = params.get('nrci_value', 0.5) # Default NRCI if not provided
        
        # Use NRCI to scale the base amplitude and overall complexity/regularity of the structure
        # Higher NRCI means more regular, coherent, or well-defined geometry
        base_amplitude_factor = 0.5 + nrci_value * 0.5 # Scales from 0.5 to 1.0 (more defined as NRCI increases)
        
        t = np.linspace(0, 2*np.pi, 100)
        vertices = []
        
        for i in range(harmonic_order):
            harmonic_freq = (i + 1) * resonance_freq # Use RGDL's primary resonance_freq for harmonics
            amplitude = base_amplitude_factor / (i + 1)  # Decreasing amplitude for higher harmonics
            
            x = amplitude * np.cos(harmonic_freq * t) * np.cos(2 * np.pi * resonance_freq * t)
            y = amplitude * np.sin(harmonic_freq * t) * np.sin(2 * np.pi * resonance_freq * t)
            z = amplitude * np.sin(2 * harmonic_freq * t) * (1.0 + (nrci_value * 0.1 * np.sin(t))) # Add NRCI influence to z-axis
            
            for j in range(len(t)):
                vertices.append([x[j], y[j], z[j]])
        
        coordinates = np.array(vertices)
        
        return GeometricPrimitive(
            primitive_type='htr_structure',
            coordinates=coordinates,
            properties={'harmonic_order': harmonic_order, 'nrci_value_influence': nrci_value},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='harmonic_toggle_resonance_nrci_influenced', # Updated method name
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    def _generate_crv_manifold(self, resonance_freq: float, params: Dict[str, Any]) -> GeometricPrimitive:
        """Generate manifold based on Core Resonance Values."""
        crv_type = params.get('crv_type', 'electromagnetic')
        manifold_dimension = params.get('manifold_dimension', 2)
        
        # Get appropriate CRV (using the internal resonance_frequencies mapping)
        crv_value = self.resonance_frequencies.get(crv_type, self.resonance_frequencies['electromagnetic'])
        
        # Generate manifold points
        if manifold_dimension == 2:
            u = np.linspace(0, 2*np.pi, 50)
            v = np.linspace(0, np.pi, 25)
            U, V = np.meshgrid(u, v)
            
            # CRV-influenced manifold
            X = np.cos(U) * np.sin(V) * (1 + 0.1 * np.sin(crv_value * U))
            Y = np.sin(U) * np.sin(V) * (1 + 0.1 * np.cos(crv_value * V))
            Z = np.cos(V) * (1 + 0.1 * np.sin(crv_value * (U + V)))
            
            vertices = []
            for i in range(X.shape[0]):
                for j in range(X.shape[1]):
                    vertices.append([X[i, j], Y[i, j], Z[i, j]])
        else:
            # 1D manifold (curve)
            t = np.linspace(0, 4*np.pi, 200)
            vertices = []
            
            for i, t_val in enumerate(t):
                x = np.cos(t_val) * (1 + 0.2 * np.sin(crv_value * t_val))
                y = np.sin(t_val) * (1 + 0.2 * np.cos(crv_val * t_val))
                z = 0.5 * np.sin(2 * t_val) * np.sin(crv_value * t_val)
                vertices.append([x, y, z])
        
        coordinates = np.array(vertices)
        
        return GeometricPrimitive(
            primitive_type='crv_manifold',
            coordinates=coordinates,
            properties={'crv_type': crv_type, 'manifold_dimension': manifold_dimension, 'crv_value': crv_value},
            resonance_frequency=resonance_freq,
            coherence_level=0.0,
            generation_method='crv_manifold_generation',
            stability_score=0.0,
            creation_timestamp=time.time()
        )
    
    # ========================================================================
    # ANALYSIS AND METRICS
    # ========================================================================
    
    def _calculate_coherence(self, primitive: GeometricPrimitive) -> float:
        """Calculate coherence level for a geometric primitive."""
        # Ensure coordinates are in the right shape for pdist and handling single points
        coords = np.array(primitive.coordinates)
        if coords.ndim == 1:
            coords = coords.reshape(1, -1) # Convert [x,y,z] to [[x,y,z]]
        
        if coords.shape[0] < 2: # For single points or empty arrays
            # For a single point, coherence is based on coordinate regularity if it has multiple dimensions,
            # otherwise it's perfectly coherent.
            if coords.shape[1] > 1: # E.g., a point [1,1,1] vs [1,2,3]
                coord_variance = np.var(coords.flatten())
                coord_mean = np.mean(np.abs(coords.flatten()))
                if coord_mean > 0:
                    coherence = 1.0 / (1.0 + coord_variance / coord_mean)
                else: # All coordinates are zero
                    coherence = 1.0
            else: # A single 1D point (scalar)
                coherence = 1.0
            return min(max(coherence, 0.0), 1.0)
        
        # Handle multi-point primitives
        try:
            distances = pdist(coords)
            if len(distances) == 0: # Should not happen if coords.shape[0] >= 2, but a safeguard
                return 1.0
            
            distance_variance = np.var(distances)
            distance_mean = np.mean(distances)
            
            # Coherence is inversely related to relative variance
            if distance_mean > 0:
                coherence = 1.0 / (1.0 + distance_variance / distance_mean)
            else: # All distances are zero (all points are identical)
                coherence = 1.0
        except Exception:
            # Fallback for any array shape or pdist issues
            coherence = 1.0
        
        return min(1.0, max(0.0, coherence))
    
    def _calculate_stability(self, primitive: GeometricPrimitive) -> float:
        """Calculate stability score for a geometric primitive."""
        coords = np.array(primitive.coordinates)
        if coords.ndim == 1:
            coords = coords.reshape(1, -1) # Ensure 2D for consistent handling
        
        num_points = coords.shape[0]
        dimension = coords.shape[1]

        # A single point or two points (a line segment) are inherently stable or 1.0
        if num_points < 3:
            return 1.0
        
        # Calculate stability based on geometric properties
        try:
            if dimension >= 3:
                # For 3D primitives, calculate convex hull volume stability
                # Need at least 4 non-coplanar points for a 3D hull
                if num_points < 4: 
                    # Fallback for insufficient points to form a 3D hull
                    # Use coordinate spread as a measure of stability
                    coord_range = np.ptp(coords, axis=0) # Peak-to-peak range for each dimension
                    stability = 1.0 / (1.0 + np.std(coord_range)) # Inverse of std of ranges
                else:
                    hull = ConvexHull(coords)
                    volume = hull.volume
                    surface_area = hull.area
                    
                    # Stability related to volume-to-surface ratio (compactness)
                    if surface_area > 0:
                        stability = volume / surface_area
                    else:
                        stability = 0.0 # Degenerate case (e.g., all points collinear/coplanar)
            elif dimension == 2:
                # For 2D primitives, calculate convex hull area stability
                # Need at least 3 non-collinear points for a 2D hull
                if num_points < 3:
                    # This case should be caught by the initial `num_points < 3` check,
                    # but kept for explicitness in this branch.
                    coord_range = np.ptp(coords, axis=0)
                    stability = 1.0 / (1.0 + np.std(coord_range))
                else:
                    # ConvexHull for 2D points, hull.volume gives the area of the polygon
                    hull = ConvexHull(coords) # Use full coordinates, ConvexHull handles 2D
                    area = hull.volume # For 2D, hull.volume is area
                    perimeter = hull.area # For 2D, hull.area is perimeter
                    
                    if perimeter > 0:
                        stability = area / perimeter
                    else:
                        stability = 0.0 # Degenerate case
            else: # 1D points or other unexpected dimensions
                coord_std = np.std(coords)
                stability = 1.0 / (1.0 + coord_std)
        except Exception as e:
            # Fallback for any other unexpected issues with ConvexHull (e.g., degenerate points)
            # print(f"DEBUG: ConvexHull failed for {primitive.primitive_type} with {num_points} points in {dimension}D: {e}")
            coord_std = np.std(coords)
            stability = 1.0 / (1.0 + coord_std)
        
        return min(1.0, max(0.0, stability))
    
    def _calculate_nrci(self, primitive: GeometricPrimitive) -> float:
        """Calculate Non-Random Coherence Index for the primitive."""
        if primitive.coordinates.size < 2: # Use .size for total number of elements
            return 1.0
        
        # Ensure coords_flat has at least 2 elements for correlation
        coords_flat = primitive.coordinates.flatten()
        if len(coords_flat) < 2:
             return 1.0 # Cannot calculate correlation for less than 2 points

        # Generate expected pattern based on resonance frequency
        # The 't' array represents a conceptual progression along the flattened coordinates.
        t = np.linspace(0, 1, len(coords_flat))
        expected_pattern = np.sin(2 * np.pi * primitive.resonance_frequency * t)
        
        # Normalize both signals to prevent bias from scale differences
        # Add a small epsilon to standard deviation to avoid division by zero if all values are identical
        epsilon = 1e-9 
        
        std_coords = np.std(coords_flat)
        if std_coords > epsilon:
            coords_normalized = (coords_flat - np.mean(coords_flat)) / std_coords
        else:
            coords_normalized = np.zeros_like(coords_flat) # If all same, normalized to zeros

        std_pattern = np.std(expected_pattern)
        if std_pattern > epsilon:
            pattern_normalized = (expected_pattern - np.mean(expected_pattern)) / std_pattern
        else:
            pattern_normalized = np.zeros_like(expected_pattern) # If all same, normalized to zeros
        
        # Calculate correlation. If one or both arrays have zero standard deviation, corrcoef returns NaN.
        # This occurs if all values are identical in an array.
        if len(coords_normalized) == len(pattern_normalized) and len(coords_normalized) > 1:
            correlation = np.corrcoef(coords_normalized, pattern_normalized)[0, 1]
            if np.isnan(correlation):
                correlation = 0.0 # Treat NaN correlation (e.g., from constant signals) as zero coherence
        else:
            correlation = 0.0 # Should not happen if len(coords_flat) >= 2 is checked, but a safeguard
        
        # Convert correlation to Non-Random Coherence Index (NRCI) (0 to 1 scale)
        nrci = (correlation + 1) / 2
        
        return min(1.0, max(0.0, nrci))
    
    def _update_metrics(self, primitive: GeometricPrimitive, generation_time: float) -> None:
        """Update performance metrics after generating a primitive."""
        self.metrics.total_primitives_generated += 1
        self.metrics.generation_time += generation_time
        
        # Update averages
        total = self.metrics.total_primitives_generated
        # Use simple moving average for clarity and performance
        self.metrics.average_coherence = (self.metrics.average_coherence * (total - 1) + 
                                        primitive.coherence_level) / total if total > 0 else 0.0
        self.metrics.average_stability = (self.metrics.average_stability * (total - 1) + 
                                        primitive.stability_score) / total if total > 0 else 0.0
        self.metrics.nrci_average = (self.metrics.nrci_average * (total - 1) + 
                                   primitive.nrci_score) / total if total > 0 else 0.0
        
        # Update resonance distribution
        realm = self._get_realm_from_frequency(primitive.resonance_frequency)
        self.metrics.resonance_distribution[realm] = self.metrics.resonance_distribution.get(realm, 0) + 1
        
        # Calculate geometric complexity (based on number of vertices)
        # Handle cases where primitive.coordinates might be empty or 1D
        num_vertices = len(primitive.coordinates) if primitive.coordinates.ndim >= 1 else 1
        complexity = num_vertices / 1000.0  # Normalize
        self.metrics.geometric_complexity = (self.metrics.geometric_complexity * (total - 1) + 
                                           complexity) / total if total > 0 else 0.0
    
    def _get_realm_from_frequency(self, frequency: float) -> str:
        """Get realm name from resonance frequency."""
        for realm, freq in self.resonance_frequencies.items():
            if abs(freq - frequency) < 1e-10: # Use a small tolerance for float comparison
                return realm
        return 'unknown'
    
    # ========================================================================
    # GEOMETRIC FIELD OPERATIONS
    # ========================================================================
    
    def create_geometric_field(self, field_name: str, primitives: List[GeometricPrimitive]) -> GeometricField:
        """
        Create a geometric field from a collection of primitives.
        
        Args:
            field_name: Name for the geometric field
            primitives: List of geometric primitives
            
        Returns:
            Created GeometricField
        """
        if not primitives:
            raise ValueError("Cannot create field with no primitives")
        
        # Calculate spatial bounds
        # Ensure all_coords has at least 2 dimensions for min/max on axis 0
        all_coords = np.vstack([p.coordinates for p in primitives])
        if all_coords.ndim == 1: # Handle case where all primitives are single points (e.g., [x,y,z])
             all_coords = all_coords.reshape(1, -1)
        
        x_bounds = (np.min(all_coords[:, 0]), np.max(all_coords[:, 0]))
        y_bounds = (np.min(all_coords[:, 1]), np.max(all_coords[:, 1]))
        z_bounds = (np.min(all_coords[:, 2]), np.max(all_coords[:, 2]))
        
        # Calculate field coherence
        coherences = [p.coherence_level for p in primitives]
        field_coherence = np.mean(coherences)
        
        # Generate resonance pattern
        resonance_pattern = np.array([p.resonance_frequency for p in primitives])
        
        # Calculate interaction matrix
        n_primitives = len(primitives)
        interaction_matrix = np.zeros((n_primitives, n_primitives))
        
        for i in range(n_primitives):
            for j in range(i + 1, n_primitives):
                # Calculate interaction strength based on resonance similarity
                freq_diff = abs(primitives[i].resonance_frequency - primitives[j].resonance_frequency)
                interaction_strength = np.exp(-freq_diff / 1000.0)  # Decay with frequency difference
                interaction_matrix[i, j] = interaction_strength
                interaction_matrix[j, i] = interaction_strength
        
        # Calculate field energy
        field_energy = np.sum([p.resonance_frequency * p.coherence_level for p in primitives])
        
        field = GeometricField(
            field_name=field_name,
            primitives=primitives,
            spatial_bounds=(x_bounds, y_bounds, z_bounds),
            field_coherence=field_coherence,
            resonance_pattern=resonance_pattern,
            interaction_matrix=interaction_matrix,
            field_energy=field_energy
        )
        
        self.geometric_fields[field_name] = field
        
        # Store in HexDictionary if available
        if self.hex_dictionary:
            field_data = {
                'field_name': field_name,
                'num_primitives': len(primitives),
                'spatial_bounds': field.spatial_bounds,
                'field_coherence': field_coherence,
                'field_energy': field_energy,
                'primitive_types': [p.primitive_type for p in primitives]
            }
            self.hex_dictionary.store(field_data, 'json', {'field_name': field_name})
        
        return field
    
    def analyze_field_interactions(self, field_name: str) -> Dict[str, Any]:
        """
        Analyze interactions within a geometric field.
        
        Args:
            field_name: Name of the field to analyze
            
        Returns:
            Dictionary with interaction analysis results
        """
        if field_name not in self.geometric_fields:
            raise ValueError(f"Field {field_name} not found")
        
        field = self.geometric_fields[field_name]
        
        # Analyze interaction matrix
        if field.interaction_matrix.size == 0: # Handle empty matrix (e.g., field with one primitive)
            return {
                'field_name': field_name,
                'average_interaction_strength': 0.0,
                'max_interaction_strength': 0.0,
                'min_interaction_strength': 0.0,
                'field_stability': 0.0,
                'strongest_interactions': [],
                'field_coherence': field.field_coherence,
                'field_energy': field.field_energy,
                'num_primitives': len(field.primitives)
            }
        
        interaction_strength = np.mean(field.interaction_matrix)
        max_interaction = np.max(field.interaction_matrix)
        min_interaction = np.min(field.interaction_matrix)
        
        # Find strongest interacting pairs
        n = field.interaction_matrix.shape[0]
        strongest_pairs = []
        
        for i in range(n):
            for j in range(i + 1, n):
                strength = field.interaction_matrix[i, j]
                strongest_pairs.append((i, j, strength))
        
        strongest_pairs.sort(key=lambda x: x[2], reverse=True)
        top_pairs = strongest_pairs[:5]  # Top 5 interactions
        
        # Calculate field stability from eigenvalues
        # Handle cases where matrix is 1x1 or cannot compute eigenvalues
        field_stability = 0.0
        if n > 0:
            try:
                eigenvalues = np.linalg.eigvals(field.interaction_matrix)
                field_stability = np.real(np.max(eigenvalues))
            except np.linalg.LinAlgError:
                field_stability = 0.0 # Cannot compute eigenvalues for singular or non-square matrix
        
        return {
            'field_name': field_name,
            'average_interaction_strength': interaction_strength,
            'max_interaction_strength': max_interaction,
            'min_interaction_strength': min_interaction,
            'field_stability': field_stability,
            'strongest_interactions': top_pairs,
            'field_coherence': field.field_coherence,
            'field_energy': field.field_energy,
            'num_primitives': len(field.primitives)
        }
    
    # ========================================================================
    # UTILITY METHODS
    # ========================================================================
    
    def get_metrics(self) -> RGDLMetrics:
        """Get current RGDL engine metrics."""
        return self.metrics
    
    def clear_cache(self) -> None:
        """Clear the geometry cache."""
        self.geometry_cache.clear()
        print("✅ RGDL geometry cache cleared")
    
    def export_primitives(self, file_path: str) -> bool:
        """
        Export all generated primitives to a file.
        
        Args:
            file_path: Path to export file
            
        Returns:
            True if successful, False otherwise
        """
        try:
            export_data = {
                'primitives': [],
                'metrics': self.metrics.__dict__,
                'geometric_fields': {}
            }
            
            # Export primitives
            for primitive in self.primitives:
                primitive_data = {
                    'primitive_type': primitive.primitive_type,
                    'coordinates': primitive.coordinates.tolist(),
                    'properties': primitive.properties,
                    'resonance_frequency': primitive.resonance_frequency,
                    'coherence_level': primitive.coherence_level,
                    'generation_method': primitive.generation_method,
                    'stability_score': primitive.stability_score,
                    'creation_timestamp': primitive.creation_timestamp,
                    'nrci_score': primitive.nrci_score
                }
                export_data['primitives'].append(primitive_data)
            
            # Export geometric fields
            for field_name, field in self.geometric_fields.items():
                field_data = {
                    'field_name': field.field_name,
                    'spatial_bounds': field.spatial_bounds,
                    'field_coherence': field.field_coherence,
                    'resonance_pattern': field.resonance_pattern.tolist(),
                    'interaction_matrix': field.interaction_matrix.tolist(),
                    'field_energy': field.field_energy,
                    'num_primitives': len(field.primitives)
                }
                export_data['geometric_fields'][field_name] = field_data
            
            with open(file_path, 'w') as f:
                json.dump(export_data, f, indent=2, default=str)
            
            print(f"✅ Exported {len(self.primitives)} primitives to {file_path}")
            return True
            
        except Exception as e:
            print(f"❌ Export failed: {e}")
            return False
    
    def generate_batch_primitives(self, batch_config: List[Dict[str, Any]]) -> List[GeometricPrimitive]:
        """
        Generate multiple primitives in batch for efficiency.
        
        Args:
            batch_config: List of configuration dictionaries for each primitive
            
        Returns:
            List of generated primitives
        """
        start_time = time.time()
        generated_primitives = []
        
        for config in batch_config:
            primitive_type = config.get('type', 'point')
            resonance_realm = config.get('realm', 'electromagnetic')
            parameters = config.get('parameters', {})
            
            try:
                primitive = self.generate_primitive(primitive_type, resonance_realm, parameters)
                generated_primitives.append(primitive)
            except Exception as e:
                print(f"⚠️  Failed to generate {primitive_type}: {e}")
                continue
        
        batch_time = time.time() - start_time
        print(f"✅ Generated {len(generated_primitives)} primitives in {batch_time:.3f}s")
        
        return generated_primitives


# ========================================================================
# UTILITY FUNCTIONS
# ========================================================================

def create_rgdl_engine(bitfield: Optional[Bitfield] = None,
                      toggle_algebra: Optional[ToggleAlgebra] = None,
                      hex_dictionary: Optional[HexDictionary] = None) -> RGDLEngine:
    """
    Create and return a new RGDL Engine instance.
    
    Args:
        bitfield: Optional Bitfield instance
        toggle_algebra: Optional ToggleAlgebra instance
        hex_dictionary: Optional HexDictionary instance
        
    Returns:
        Initialized RGDLEngine instance
    """
    return RGDLEngine(bitfield, toggle_algebra, hex_dictionary)


def benchmark_rgdl_engine(engine: RGDLEngine, num_primitives: int = 100) -> Dict[str, float]:
    """
    Benchmark RGDL Engine performance.
    
    Args:
        engine: RGDLEngine instance to benchmark
        num_primitives: Number of primitives to generate
        
    Returns:
        Dictionary with benchmark results
    """
    start_time = time.time()
    
    # Generate various primitive types
    primitive_types = ['point', 'line', 'triangle', 'sphere', 'cube']
    realms = ['quantum', 'electromagnetic', 'gravitational', 'biological']
    
    generation_times = []
    
    for i in range(num_primitives):
        primitive_type = primitive_types[i % len(primitive_types)]
        realm = realms[i % len(realms)]
        
        gen_start = time.time()
        engine.generate_primitive(primitive_type, realm)
        generation_times.append(time.time() - gen_start)
    
    total_time = time.time() - start_time
    
    return {
        'total_time': total_time,
        'average_generation_time': np.mean(generation_times),
        'primitives_per_second': num_primitives / total_time,
        'total_primitives_generated': engine.metrics.total_primitives_generated,
        'average_coherence': engine.metrics.average_coherence,
        'average_stability': engine.metrics.average_stability,
        'average_nrci': engine.metrics.nrci_average
    }


if __name__ == "__main__":
    # Test the RGDL Engine
    print("🧪 Testing RGDL Engine v3.1...")
    
    # To run this test independently, ensure UBPConstants and HexDictionary are available
    # For a full system test, this should be run via _run_ubp_test.py
    
    # Minimal mock for UBPConstants if not run within full UBP environment
    class MockUBPConstants:
        CRV_QUANTUM_BASE = 1.6180339887  # Golden ratio, example value
        CRV_ELECTROMAGNETIC_BASE = 2.7182818284 # Euler's number, example value
        CRV_GRAVITATIONAL_BASE = 3.1415926535 # Pi, example value
        CRV_BIOLOGICAL_BASE = 0.5772156649 # Euler-Mascheroni constant, example value
        CRV_COSMOLOGICAL_BASE = 1.0
        CRV_NUCLEAR_BASE = 137.035999139 # Fine-structure constant inverse, example value
        CRV_OPTICAL_BASE = 299792458 # Speed of light, example value

    # Override UBPConstants for standalone testing if needed
    try:
        # Changed core to system_constants as constants are now defined there
        from system_constants import UBPConstants 
    except ImportError:
        UBPConstants = MockUBPConstants()
        print("⚠️ Using MockUBPConstants for standalone test.")

    # Mock HexDictionary for standalone testing
    class MockHexDictionary:
        def store(self, data, data_format, metadata):
            # print(f"MockHexDictionary storing: {metadata['field_name'] if 'field_name' in metadata else metadata['primitive_type']}")
            pass
        def retrieve(self, query):
            return None

    engine = create_rgdl_engine(hex_dictionary=MockHexDictionary())
    
    # Test basic primitive generation
    point = engine.generate_primitive('point', 'quantum')
    line = engine.generate_primitive('line', 'electromagnetic')
    sphere = engine.generate_primitive('sphere', 'gravitational')
    
    print(f"Generated point: coherence={point.coherence_level:.3f}, NRCI={point.nrci_score:.3f}")
    print(f"Generated line: coherence={line.coherence_level:.3f}, NRCI={line.nrci_score:.3f}")
    print(f"Generated sphere: coherence={sphere.coherence_level:.3f}, NRCI={sphere.nrci_score:.3f}")
    
    # Test geometric field creation
    field = engine.create_geometric_field('test_field', [point, line, sphere])
    field_analysis = engine.analyze_field_interactions('test_field')
    
    print(f"Field coherence: {field_analysis['field_coherence']:.3f}")
    print(f"Field energy: {field_analysis['field_energy']:.3f}")
    
    # Test metrics
    metrics = engine.get_metrics()
    print(f"Total primitives: {metrics.total_primitives_generated}")
    print(f"Average coherence: {metrics.average_coherence:.3f}")
    print(f"Average NRCI: {metrics.nrci_average:.3f}")
    
    print("✅ RGDL Engine v3.1 test completed successfully!")