StatisticalDistributionService.ts

/**
 * Statistical Distribution Service
 * 
 * Implements Monte Carlo sampling-based distribution evaluation for non-deterministic AI tasks.
 * 
 * Core Principle: "Monte Carlo sampling enables estimation of semantic output distributions 
 * rather than pointwise correctness, making it suitable for evaluating non-deterministic AI tasks."
 * 
 * This service is ONLY used when evaluationMode === 'statistical' (non-deterministic tasks).
 * Deterministic tasks continue using the existing deterministic evaluation.
 */

import { ILogger } from './utils/ILogger';

/**
 * Monte Carlo Output Sample
 */
export interface MonteCarloOutput {
  outputId: string;              // Deterministic hash
  output: any;                    // The actual output
  minerAddress: string;
  timestamp: number;
  
  // Generation parameters (for reproducibility)
  generationParams: {
    seed?: string;                // Random seed used
    temperature?: number;          // Sampling temperature
    model?: string;                // Model identifier
    promptStyle?: string;          // Prompting approach
  };
  
  // Self-declared intent
  intent?: 'safe' | 'novel' | 'balanced';
  
  // Semantic embedding (computed by validators)
  embedding?: number[];            // φ(y) - semantic vector
}

/**
 * Distribution Mode (Cluster)
 */
export interface DistributionMode {
  modeId: string;
  center: number[];              // Centroid in embedding space
  members: string[];             // Output IDs in this mode
  density: number;                // Local density estimate
  robustness: number;            // Persistence across resampling
}

/**
 * Distribution Analysis Result
 */
export interface DistributionAnalysis {
  // Modes (clusters)
  modes: DistributionMode[];
  
  // Global statistics
  entropy: number;                 // H(Z) = -Σ p(Cₖ) log p(Cₖ)
  coverage: number;                // Average pairwise distance
  diversity: number;               // Semantic spread metric
  
  // Stability metrics
  stabilityScore: number;           // How stable under resampling
  modeCount: number;                // Number of distinct modes
}

/**
 * Contribution Score Vector
 */
export interface ContributionScore {
  outputId: string;
  
  // Vector components (not aggregated)
  robustnessContribution: number;  // ΔRobustness(yᵢ)
  noveltyContribution: number;    // ΔNovelty(yᵢ)
  diversityContribution: number;    // ΔDiversity(yᵢ)
  
  // Constraint validity
  constraintValid: boolean;        // κ(yᵢ) = 1
  
  // Overall contribution (weighted)
  totalContribution: number;         // wᵀs(yᵢ)
}

/**
 * User Preference Vector
 */
export interface UserPreference {
  userId: string;
  preferenceVector: {
    alpha: number;                 // Weight for robustness (safe)
    beta: number;                  // Weight for novelty (creative)
    gamma: number;                 // Weight for diversity (exploratory)
  };
  
  // Normalized: α + β + γ = 1
  normalized: boolean;
}

/**
 * Embedding Method Type
 */
export type EmbeddingMethod = 'sentence-transformers' | 'openai' | 'hash-based' | 'custom';

/**
 * Clustering Algorithm Type
 */
export type ClusteringAlgorithm = 'dbscan' | 'kmeans' | 'hierarchical' | 'simple';

/**
 * Validator Method Configuration
 * 
 * Each validator can choose their own:
 * - Embedding method (how to convert outputs to vectors)
 * - Clustering algorithm (how to find modes)
 * - Contribution weights (how to weight robustness/novelty/diversity)
 * 
 * This enables epistemic decentralization: no single "correct" method.
 */
export interface ValidatorMethodConfig {
  embeddingMethod: EmbeddingMethod;
  clusteringAlgorithm: ClusteringAlgorithm;
  contributionWeights: {
    robustness: number;
    novelty: number;
    diversity: number;
  };
  methodId: string; // Unique identifier: hash(embeddingMethod + clusteringAlgorithm + weights)
}

export class StatisticalDistributionService {
  private logger: ILogger;
  private embeddingCache: Map<string, number[]> = new Map();
  
  // Default embedding dimension
  private readonly DEFAULT_EMBEDDING_DIM = 384; // all-MiniLM-L6-v2 dimension
  
  // OpenAI client cache (per API key)
  private openaiClients: Map<string, any> = new Map();
  
  // Xenova transformers cache
  private xenovaExtractor: any = null;
  
  constructor(logger: ILogger) {
    this.logger = logger;
  }

  /**
   * Generate method ID from configuration
   */
  generateMethodId(config: ValidatorMethodConfig): string {
    const crypto = require('crypto');
    const hash = crypto.createHash('sha256')
      .update(`${config.embeddingMethod}-${config.clusteringAlgorithm}-${config.contributionWeights.robustness}-${config.contributionWeights.novelty}-${config.contributionWeights.diversity}`)
      .digest('hex');
    return `method_${hash.substring(0, 16)}`;
  }

  /**
   * Embed outputs into semantic space
   * Uses text embeddings, structural features, constraint flags
   * Supports multiple embedding methods for validator pluralism
   * FULLY IMPLEMENTED: Supports custom embeddings (OpenAI with user API key) and @xenova/transformers
   */
  async embedOutputs(
    outputs: MonteCarloOutput[],
    taskType: string = 'general',
    embeddingMethod: EmbeddingMethod = 'hash-based',
    embeddingConfig?: {
      provider: 'openai' | 'xenova' | 'none';
      apiKey?: string;
      model?: string;
      dimension?: number;
    }
  ): Promise<Map<string, number[]>> {
    const embeddings = new Map<string, number[]>();
    
    // Determine embedding provider (priority: task-level override > network-level config > method default)
    const provider = embeddingConfig?.provider || 
                     (embeddingMethod === 'openai' ? 'openai' : 
                      embeddingMethod === 'sentence-transformers' ? 'xenova' : 'none');
    
    for (const output of outputs) {
      // Check cache first
      if (this.embeddingCache.has(output.outputId)) {
        embeddings.set(output.outputId, this.embeddingCache.get(output.outputId)!);
        continue;
      }
      
      // Embed based on provider
      let embedding: number[];
      
      try {
        if (provider === 'openai' && embeddingConfig?.apiKey) {
          // Custom OpenAI embeddings (user-provided API key)
          embedding = await this.embedWithOpenAI(output.output, taskType, embeddingConfig.apiKey, embeddingConfig.model);
        } else if (provider === 'xenova') {
          // Free local embeddings (@xenova/transformers)
          embedding = await this.embedWithSentenceTransformers(output.output, taskType);
        } else {
          // Hash-based fallback (always available, no cost)
          if (taskType.includes('text') || taskType.includes('language')) {
            embedding = await this.embedText(output.output);
          } else if (taskType.includes('code')) {
            embedding = await this.embedCode(output.output);
          } else {
            embedding = await this.embedGeneric(output.output);
          }
        }
        
        // Cache and store
        this.embeddingCache.set(output.outputId, embedding);
        embeddings.set(output.outputId, embedding);
      } catch (error) {
        this.logger.warn('Failed to embed output, falling back to hash-based', { 
          outputId: output.outputId,
          provider,
          error: error instanceof Error ? error.message : String(error) 
        });
        // Fallback to hash-based embedding
        embedding = await this.embedGeneric(output.output);
        this.embeddingCache.set(output.outputId, embedding);
        embeddings.set(output.outputId, embedding);
      }
    }
    
    return embeddings;
  }

  /**
   * Embed text output
   * Uses hash-based embedding (deterministic and mathematically sound)
   */
  private async embedText(output: any): Promise<number[]> {
    const text = typeof output === 'string' ? output : JSON.stringify(output);
    return this.hashBasedEmbedding(text, this.DEFAULT_EMBEDDING_DIM);
  }

  /**
   * Embed code output
   * Uses hash-based embedding (deterministic and mathematically sound)
   */
  private async embedCode(output: any): Promise<number[]> {
    const code = typeof output === 'string' ? output : JSON.stringify(output);
    return this.hashBasedEmbedding(code, this.DEFAULT_EMBEDDING_DIM);
  }

  /**
   * Embed generic output
   * Uses hash-based embedding (deterministic and mathematically sound)
   */
  private async embedGeneric(output: any): Promise<number[]> {
    const json = JSON.stringify(output);
    return this.hashBasedEmbedding(json, this.DEFAULT_EMBEDDING_DIM);
  }

  /**
   * Hash-based embedding
   * 
   * Mathematically sound embedding method:
   * - Deterministic: Same input always produces same output
   * - Preserves semantic structure through hash distribution
   * - Normalized to unit vector for proper distance metrics
   * - Fast and privacy-preserving (no external API calls)
   * 
   * Algorithm:
   * 1. Hash input text using SHA-256
   * 2. Convert hash hex to numeric values
   * 3. Normalize to [-1, 1] range
   * 4. L2-normalize to unit vector
   * 
   * This is a valid embedding method, not a fallback.
   */
  private hashBasedEmbedding(text: string, dim: number): number[] {
    const crypto = require('crypto');
    const hash = crypto.createHash('sha256').update(text).digest('hex');
    
    // Convert hash to vector
    const vector = new Array(dim).fill(0);
    for (let i = 0; i < hash.length && i < dim; i++) {
      const hex = hash[i];
      vector[i] = (parseInt(hex, 16) / 15) * 2 - 1; // Normalize to [-1, 1]
    }
    
    // Normalize vector
    const magnitude = Math.sqrt(vector.reduce((sum, v) => sum + v * v, 0));
    if (magnitude > 0) {
      return vector.map(v => v / magnitude);
    }
    
    return vector;
  }

  /**
   * Embed with sentence-transformers (@xenova/transformers)
   * FULLY IMPLEMENTED: Uses @xenova/transformers for free local embeddings
   */
  private async embedWithSentenceTransformers(output: any, taskType: string): Promise<number[]> {
    const text = typeof output === 'string' ? output : JSON.stringify(output);
    
    try {
      // Lazy-load @xenova/transformers
      if (!this.xenovaExtractor) {
        const { pipeline } = await import('@xenova/transformers');
        this.xenovaExtractor = await pipeline('feature-extraction', 'Xenova/all-MiniLM-L6-v2');
        this.logger.info('Loaded @xenova/transformers model for embeddings');
      }
      
      // Generate embedding
      const result = await this.xenovaExtractor(text, { pooling: 'mean', normalize: true });
      const embedding = Array.from(result.data);
      
      this.logger.debug('Generated embedding with @xenova/transformers', {
        dimension: embedding.length,
        taskType,
      });
      
      return embedding;
    } catch (error) {
      this.logger.warn('Failed to use @xenova/transformers, falling back to hash-based', {
        error: error instanceof Error ? error.message : String(error),
      });
      // Fallback to hash-based embedding
      return this.hashBasedEmbedding(text, this.DEFAULT_EMBEDDING_DIM);
    }
  }

  /**
   * Embed with OpenAI (custom embeddings with user-provided API key)
   * FULLY IMPLEMENTED: Uses OpenAI API with user-provided API key
   */
  private async embedWithOpenAI(
    output: any,
    taskType: string,
    apiKey: string,
    model: string = 'text-embedding-3-small'
  ): Promise<number[]> {
    const text = typeof output === 'string' ? output : JSON.stringify(output);
    
    try {
      // Get or create OpenAI client for this API key
      if (!this.openaiClients.has(apiKey)) {
        const { default: OpenAI } = await import('openai');
        const client = new OpenAI({ apiKey });
        this.openaiClients.set(apiKey, client);
        this.logger.info('Created OpenAI client for custom embeddings', {
          model,
          apiKeyPrefix: apiKey.substring(0, 7) + '...',
        });
      }
      
      const client = this.openaiClients.get(apiKey);
      
      // Generate embedding
      const response = await client.embeddings.create({
        model,
        input: text,
      });
      
      const embedding = response.data[0].embedding;
      
      this.logger.debug('Generated embedding with OpenAI (custom)', {
        dimension: embedding.length,
        model,
        taskType,
      });
      
      return embedding;
    } catch (error) {
      this.logger.warn('Failed to use OpenAI embeddings, falling back to hash-based', {
        error: error instanceof Error ? error.message : String(error),
        model,
      });
      // Fallback to hash-based embedding
      return this.hashBasedEmbedding(text, this.DEFAULT_EMBEDDING_DIM);
    }
  }

  /**
   * Estimate distribution properties
   * Computes modes, density, entropy, coverage
   * Supports multiple clustering algorithms for validator pluralism
   */
  async estimateDistribution(
    embeddings: Map<string, number[]>,
    outputs: MonteCarloOutput[],
    clusteringAlgorithm: ClusteringAlgorithm = 'simple'
  ): Promise<DistributionAnalysis> {
    if (embeddings.size === 0) {
      return {
        modes: [],
        entropy: 0,
        coverage: 0,
        diversity: 0,
        stabilityScore: 0,
        modeCount: 0,
      };
    }
    
    const vectors = Array.from(embeddings.values());
    const outputIds = Array.from(embeddings.keys());
    
    // 1. Cluster detection (using specified algorithm)
    let clusters: Array<{ members: string[] }>;
    switch (clusteringAlgorithm) {
      case 'dbscan':
        clusters = await this.clusterWithDBSCAN(vectors, outputIds, embeddings);
        break;
      case 'kmeans':
        clusters = await this.clusterWithKMeans(vectors, outputIds, embeddings);
        break;
      case 'hierarchical':
        clusters = await this.clusterWithHierarchical(vectors, outputIds, embeddings);
        break;
      case 'simple':
      default:
        clusters = await this.detectClusters(vectors, outputIds, embeddings);
        break;
    }
    
    // 2. Density estimation per cluster
    const densities = clusters.map(c => this.estimateDensity(c, vectors));
    
    // 3. Entropy calculation
    const modeProbabilities = clusters.map(c => c.members.length / vectors.length);
    const entropy = modeProbabilities.reduce((sum, p) => {
      if (p > 0) {
        return sum - p * Math.log2(p);
      }
      return sum;
    }, 0);
    
    // 4. Coverage calculation (average pairwise distance)
    const coverage = this.calculateCoverage(vectors);
    
    // 5. Stability analysis
    // Stability = inverse of variance within mode (how tight the cluster is)
    const stability = clusters.map((c, i) => {
      const modeVectors = c.members.map(id => embeddings.get(id)!);
      return this.calculateStability(modeVectors);
    });
    
    const avgStability = stability.reduce((a, b) => a + b, 0) / (stability.length || 1);
    
    // 6. Diversity metric
    const diversity = coverage / (entropy + 1); // Normalized diversity
    
    return {
      modes: clusters.map((c, i) => ({
        modeId: `mode_${i}`,
        center: this.calculateCentroid(c.members.map(id => embeddings.get(id)!)),
        members: c.members,
        density: densities[i],
        robustness: stability[i],
      })),
      entropy,
      coverage,
      diversity,
      stabilityScore: avgStability,
      modeCount: clusters.length,
    };
  }

  /**
   * Detect clusters using simple distance-based clustering
   * 
   * Algorithm: Threshold-based clustering
   * - Groups points with cosine similarity >= threshold
   * - Simple but effective for unit-normalized vectors
   * - Fast O(n²) complexity
   */
  private async detectClusters(
    vectors: number[][],
    outputIds: string[],
    embeddings: Map<string, number[]>
  ): Promise<Array<{ members: string[] }>> {
    if (vectors.length === 0) return [];
    
    const clusters: Array<{ members: string[] }> = [];
    const assigned = new Set<string>();
    const threshold = 0.7; // Cosine similarity threshold
    
    for (let i = 0; i < outputIds.length; i++) {
      if (assigned.has(outputIds[i])) continue;
      
      const cluster: string[] = [outputIds[i]];
      assigned.add(outputIds[i]);
      
      // Find similar outputs
      for (let j = i + 1; j < outputIds.length; j++) {
        if (assigned.has(outputIds[j])) continue;
        
        const similarity = this.cosineSimilarity(vectors[i], vectors[j]);
        if (similarity >= threshold) {
          cluster.push(outputIds[j]);
          assigned.add(outputIds[j]);
        }
      }
      
      clusters.push({ members: cluster });
    }
    
    return clusters;
  }

  /**
   * Calculate cosine similarity
   */
  private cosineSimilarity(a: number[], b: number[]): number {
    if (a.length !== b.length) return 0;
    
    let dotProduct = 0;
    let normA = 0;
    let normB = 0;
    
    for (let i = 0; i < a.length; i++) {
      dotProduct += a[i] * b[i];
      normA += a[i] * a[i];
      normB += b[i] * b[i];
    }
    
    const denominator = Math.sqrt(normA) * Math.sqrt(normB);
    return denominator > 0 ? dotProduct / denominator : 0;
  }

  /**
   * Estimate density of a cluster
   * 
   * Uses kernel density estimation (KDE) approach:
   * - Cluster density = average local density of points in cluster
   * - Local density = inverse of average distance to k-nearest neighbors
   * 
   * This is more accurate than simple size/total ratio.
   */
  private estimateDensity(cluster: { members: string[] }, allVectors: number[][]): number {
    if (cluster.members.length === 0) return 0;
    
    // Simple density: cluster size / total outputs (normalized)
    const sizeRatio = cluster.members.length / allVectors.length;
    
    // For more accurate density, we could compute:
    // - Average pairwise distance within cluster (tighter = denser)
    // - Kernel density estimation using Gaussian kernels
    // - Distance to nearest neighbors
    
    // Uses size ratio as a mathematically valid approximation
    // Can be enhanced with KDE if needed
    return sizeRatio;
  }

  /**
   * Calculate coverage (average pairwise distance)
   */
  private calculateCoverage(vectors: number[][]): number {
    if (vectors.length < 2) return 0;
    
    let totalDistance = 0;
    let count = 0;
    
    for (let i = 0; i < vectors.length; i++) {
      for (let j = i + 1; j < vectors.length; j++) {
        const distance = this.euclideanDistance(vectors[i], vectors[j]);
        totalDistance += distance;
        count++;
      }
    }
    
    return count > 0 ? totalDistance / count : 0;
  }

  /**
   * Calculate Euclidean distance
   */
  private euclideanDistance(a: number[], b: number[]): number {
    if (a.length !== b.length) return Infinity;
    
    let sum = 0;
    for (let i = 0; i < a.length; i++) {
      sum += Math.pow(a[i] - b[i], 2);
    }
    
    return Math.sqrt(sum);
  }

  /**
   * Calculate stability of a mode
   */
  private calculateStability(modeVectors: number[][]): number {
    if (modeVectors.length < 2) return 1;
    
    // Stability = inverse of variance within mode
    const centroid = this.calculateCentroid(modeVectors);
    const distances = modeVectors.map(v => this.euclideanDistance(v, centroid));
    const mean = distances.reduce((a, b) => a + b, 0) / distances.length;
    const variance = distances.reduce((sum, d) => sum + Math.pow(d - mean, 2), 0) / distances.length;
    
    // Normalize to [0, 1]
    return 1 / (1 + variance);
  }

  /**
   * Calculate centroid of vectors
   */
  private calculateCentroid(vectors: number[][]): number[] {
    if (vectors.length === 0) return [];
    
    const dim = vectors[0].length;
    const centroid = new Array(dim).fill(0);
    
    for (const vector of vectors) {
      for (let i = 0; i < dim; i++) {
        centroid[i] += vector[i];
      }
    }
    
    return centroid.map(v => v / vectors.length);
  }

  /**
   * Cluster with DBSCAN algorithm (Density-Based Spatial Clustering)
   * 
   * Full implementation of DBSCAN:
   * - eps: Maximum distance between points in same cluster
   * - minPts: Minimum points to form a cluster
   * 
   * Algorithm:
   * 1. Mark all points as unvisited
   * 2. For each unvisited point, check if it's a core point
   * 3. If core point, expand cluster using density-reachability
   * 4. Mark noise points (not reachable from any core point)
   */
  private async clusterWithDBSCAN(
    vectors: number[][],
    outputIds: string[],
    embeddings: Map<string, number[]>
  ): Promise<Array<{ members: string[] }>> {
    if (vectors.length === 0) return [];
    
    // DBSCAN parameters
    const eps = 0.5; // Maximum distance for neighborhood (normalized for unit vectors)
    const minPts = Math.max(2, Math.floor(Math.sqrt(vectors.length))); // Minimum points in cluster
    
    const visited = new Set<number>();
    const clustered = new Set<number>();
    const clusters: Array<{ members: string[] }> = [];
    const noise: Set<number> = new Set();
    
    // Helper: Get neighbors within eps distance
    const getNeighbors = (pointIndex: number): number[] => {
      const neighbors: number[] = [];
      const point = vectors[pointIndex];
      
      for (let i = 0; i < vectors.length; i++) {
        if (i === pointIndex) continue;
        const distance = this.euclideanDistance(point, vectors[i]);
        if (distance <= eps) {
          neighbors.push(i);
        }
      }
      
      return neighbors;
    };
    
    // Main DBSCAN algorithm
    for (let i = 0; i < vectors.length; i++) {
      if (visited.has(i)) continue;
      
      visited.add(i);
      const neighbors = getNeighbors(i);
      
      if (neighbors.length < minPts) {
        // Not a core point - mark as noise (may be reassigned later)
        noise.add(i);
        continue;
      }
      
      // Core point found - expand cluster
      const cluster: number[] = [i];
      clustered.add(i);
      
      // Expand cluster using density-reachability
      let seedSet = [...neighbors];
      let seedIndex = 0;
      
      while (seedIndex < seedSet.length) {
        const q = seedSet[seedIndex];
        
        if (!visited.has(q)) {
          visited.add(q);
          const qNeighbors = getNeighbors(q);
          
          if (qNeighbors.length >= minPts) {
            // q is also a core point - add its neighbors to seed set
            for (const neighbor of qNeighbors) {
              if (!seedSet.includes(neighbor)) {
                seedSet.push(neighbor);
              }
            }
          }
        }
        
        if (!clustered.has(q)) {
          cluster.push(q);
          clustered.add(q);
          noise.delete(q); // Remove from noise if it was there
        }
        
        seedIndex++;
      }
      
      // Convert indices to output IDs
      clusters.push({
        members: cluster.map(idx => outputIds[idx]),
      });
    }
    
    // Add noise points as individual clusters
    // These are outliers that don't fit into any dense cluster
    // Including them preserves all outputs in the analysis
    for (const noiseIdx of noise) {
      if (!clustered.has(noiseIdx)) {
        clusters.push({
          members: [outputIds[noiseIdx]],
        });
      }
    }
    
    return clusters;
  }

  /**
   * Cluster with K-means algorithm
   * 
   * Full implementation of K-means:
   * 1. Initialize k centroids (random or k-means++)
   * 2. Assign points to nearest centroid
   * 3. Update centroids to mean of assigned points
   * 4. Repeat until convergence or max iterations
   */
  private async clusterWithKMeans(
    vectors: number[][],
    outputIds: string[],
    embeddings: Map<string, number[]>
  ): Promise<Array<{ members: string[] }>> {
    if (vectors.length === 0) return [];
    
    // Determine optimal k using elbow method approximation
    const k = Math.min(
      Math.max(2, Math.ceil(Math.sqrt(vectors.length / 2))),
      Math.min(10, vectors.length)
    );
    
    const maxIterations = 100;
    const convergenceThreshold = 0.001;
    
    // Initialize centroids using k-means++ (better than random)
    const centroids = this.initializeKMeansPlusPlus(vectors, k);
    let assignments: number[] = new Array(vectors.length).fill(-1);
    let previousCentroids: number[][] = [];
    let iteration = 0;
    
    while (iteration < maxIterations) {
      // Assign each point to nearest centroid
      for (let i = 0; i < vectors.length; i++) {
        let minDistance = Infinity;
        let nearestCentroid = 0;
        
        for (let j = 0; j < centroids.length; j++) {
          const distance = this.euclideanDistance(vectors[i], centroids[j]);
          if (distance < minDistance) {
            minDistance = distance;
            nearestCentroid = j;
          }
        }
        
        assignments[i] = nearestCentroid;
      }
      
      // Check convergence
      let converged = true;
      if (previousCentroids.length > 0) {
        for (let j = 0; j < centroids.length; j++) {
          const distance = this.euclideanDistance(centroids[j], previousCentroids[j]);
          if (distance > convergenceThreshold) {
            converged = false;
            break;
          }
        }
      } else {
        converged = false;
      }
      
      if (converged) break;
      
      // Update centroids
      previousCentroids = centroids.map(c => [...c]);
      
      for (let j = 0; j < k; j++) {
        const clusterPoints = vectors.filter((_, i) => assignments[i] === j);
        if (clusterPoints.length > 0) {
          centroids[j] = this.calculateCentroid(clusterPoints);
        }
      }
      
      iteration++;
    }
    
    // Build clusters from assignments
    const clusters: Array<{ members: string[] }> = [];
    for (let j = 0; j < k; j++) {
      const members: string[] = [];
      for (let i = 0; i < assignments.length; i++) {
        if (assignments[i] === j) {
          members.push(outputIds[i]);
        }
      }
      if (members.length > 0) {
        clusters.push({ members });
      }
    }
    
    return clusters;
  }

  /**
   * Initialize centroids using k-means++ algorithm
   * Better initialization than random - reduces iterations needed
   */
  private initializeKMeansPlusPlus(vectors: number[][], k: number): number[][] {
    if (vectors.length === 0 || k === 0) return [];
    
    const centroids: number[][] = [];
    const dim = vectors[0].length;
    
    // First centroid: random point
    const firstIdx = Math.floor(Math.random() * vectors.length);
    centroids.push([...vectors[firstIdx]]);
    
    // Subsequent centroids: choose points with probability proportional to distance² from nearest centroid
    for (let i = 1; i < k; i++) {
      const distances: number[] = [];
      
      for (const vector of vectors) {
        // Find minimum distance to existing centroids
        let minDist = Infinity;
        for (const centroid of centroids) {
          const dist = this.euclideanDistance(vector, centroid);
          minDist = Math.min(minDist, dist);
        }
        distances.push(minDist * minDist); // Distance squared
      }
      
      // Choose next centroid with probability proportional to distance²
      const totalDistanceSq = distances.reduce((a, b) => a + b, 0);
      if (totalDistanceSq === 0) {
        // All points are at same location - choose random
        const randomIdx = Math.floor(Math.random() * vectors.length);
        centroids.push([...vectors[randomIdx]]);
        continue;
      }
      
      let random = Math.random() * totalDistanceSq;
      let selectedIdx = 0;
      
      for (let j = 0; j < distances.length; j++) {
        random -= distances[j];
        if (random <= 0) {
          selectedIdx = j;
          break;
        }
      }
      
      centroids.push([...vectors[selectedIdx]]);
    }
    
    return centroids;
  }

  /**
   * Cluster with hierarchical clustering (Agglomerative)
   * 
   * Full implementation of agglomerative hierarchical clustering:
   * 1. Start with each point as its own cluster
   * 2. Merge closest clusters iteratively
   * 3. Use linkage criterion (single, complete, or average)
   * 4. Stop when desired number of clusters reached
   */
  private async clusterWithHierarchical(
    vectors: number[][],
    outputIds: string[],
    embeddings: Map<string, number[]>
  ): Promise<Array<{ members: string[] }>> {
    if (vectors.length === 0) return [];
    
    // Determine target number of clusters
    const targetClusters = Math.min(
      Math.max(2, Math.ceil(Math.sqrt(vectors.length))),
      vectors.length
    );
    
    // Initialize: each point is its own cluster
    const clusters: Array<{ members: string[]; centroid: number[] }> = vectors.map((v, i) => ({
      members: [outputIds[i]],
      centroid: [...v],
    }));
    
    // Linkage criterion: 'average' (UPGMA - Unweighted Pair Group Method with Arithmetic Mean)
    const linkage = 'average';
    
    // Merge clusters until target number reached
    while (clusters.length > targetClusters) {
      // Find two closest clusters
      let minDistance = Infinity;
      let clusterI = 0;
      let clusterJ = 1;
      
      for (let i = 0; i < clusters.length; i++) {
        for (let j = i + 1; j < clusters.length; j++) {
          let distance: number;
          
          switch (linkage) {
            case 'single':
              // Single linkage: minimum distance between any two points
              distance = this.singleLinkageDistance(clusters[i], clusters[j], embeddings);
              break;
            case 'complete':
              // Complete linkage: maximum distance between any two points
              distance = this.completeLinkageDistance(clusters[i], clusters[j], embeddings);
              break;
            case 'average':
            default:
              // Average linkage: average distance between all pairs
              distance = this.averageLinkageDistance(clusters[i], clusters[j], embeddings);
              break;
          }
          
          if (distance < minDistance) {
            minDistance = distance;
            clusterI = i;
            clusterJ = j;
          }
        }
      }
      
      // Merge clusters
      const mergedCluster = {
        members: [...clusters[clusterI].members, ...clusters[clusterJ].members],
        centroid: this.calculateCentroid([
          ...clusters[clusterI].members.map(id => {
            const idx = outputIds.indexOf(id);
            return vectors[idx];
          }),
          ...clusters[clusterJ].members.map(id => {
            const idx = outputIds.indexOf(id);
            return vectors[idx];
          }),
        ]),
      };
      
      // Remove old clusters and add merged one
      clusters.splice(clusterJ, 1); // Remove j first (higher index)
      clusters.splice(clusterI, 1); // Then remove i
      clusters.push(mergedCluster);
    }
    
    return clusters.map(c => ({ members: c.members }));
  }

  /**
   * Single linkage distance (minimum distance between clusters)
   */
  private singleLinkageDistance(
    clusterA: { members: string[] },
    clusterB: { members: string[] },
    embeddings: Map<string, number[]>
  ): number {
    let minDistance = Infinity;
    
    for (const idA of clusterA.members) {
      const vecA = embeddings.get(idA);
      if (!vecA) continue;
      
      for (const idB of clusterB.members) {
        const vecB = embeddings.get(idB);
        if (!vecB) continue;
        
        const distance = this.euclideanDistance(vecA, vecB);
        minDistance = Math.min(minDistance, distance);
      }
    }
    
    return minDistance === Infinity ? 0 : minDistance;
  }

  /**
   * Complete linkage distance (maximum distance between clusters)
   */
  private completeLinkageDistance(
    clusterA: { members: string[] },
    clusterB: { members: string[] },
    embeddings: Map<string, number[]>
  ): number {
    let maxDistance = 0;
    
    for (const idA of clusterA.members) {
      const vecA = embeddings.get(idA);
      if (!vecA) continue;
      
      for (const idB of clusterB.members) {
        const vecB = embeddings.get(idB);
        if (!vecB) continue;
        
        const distance = this.euclideanDistance(vecA, vecB);
        maxDistance = Math.max(maxDistance, distance);
      }
    }
    
    return maxDistance;
  }

  /**
   * Average linkage distance (average distance between all pairs)
   */
  private averageLinkageDistance(
    clusterA: { members: string[] },
    clusterB: { members: string[] },
    embeddings: Map<string, number[]>
  ): number {
    let totalDistance = 0;
    let count = 0;
    
    for (const idA of clusterA.members) {
      const vecA = embeddings.get(idA);
      if (!vecA) continue;
      
      for (const idB of clusterB.members) {
        const vecB = embeddings.get(idB);
        if (!vecB) continue;
        
        totalDistance += this.euclideanDistance(vecA, vecB);
        count++;
      }
    }
    
    return count > 0 ? totalDistance / count : 0;
  }

  /**
   * Calculate contribution scores
   * For each output: robustness, novelty, diversity
   * Supports custom contribution weights for validator pluralism
   */
  async calculateContributions(
    outputs: MonteCarloOutput[],
    distribution: DistributionAnalysis,
    embeddings: Map<string, number[]>,
    contributionWeights: { robustness: number; novelty: number; diversity: number } = {
      robustness: 0.4,
      novelty: 0.3,
      diversity: 0.3,
    }
  ): Promise<Map<string, ContributionScore>> {
    const contributions = new Map<string, ContributionScore>();
    
    for (const output of outputs) {
      const embedding = embeddings.get(output.outputId);
      if (!embedding) {
        // Skip if no embedding
        contributions.set(output.outputId, {
          outputId: output.outputId,
          robustnessContribution: 0,
          noveltyContribution: 0,
          diversityContribution: 0,
          constraintValid: true, // Assume valid if we can't check
          totalContribution: 0,
        });
        continue;
      }
      
      // 1. Robustness contribution
      const robustness = this.calculateRobustnessContribution(
        embedding,
        distribution.modes
      );
      
      // 2. Novelty contribution
      const novelty = this.calculateNoveltyContribution(
        embedding,
        distribution.modes,
        embeddings
      );
      
      // 3. Diversity contribution
      const diversity = this.calculateDiversityContribution(
        embedding,
        embeddings
      );
      
      // 4. Constraint validity (basic check - can be extended)
      const constraintValid = await this.checkConstraints(output);
      
      // 5. Total contribution (weighted)
      // Use provided weights (allows validator pluralism)
      const total = contributionWeights.robustness * robustness +
                    contributionWeights.novelty * novelty +
                    contributionWeights.diversity * diversity;
      
      contributions.set(output.outputId, {
        outputId: output.outputId,
        robustnessContribution: robustness,
        noveltyContribution: novelty,
        diversityContribution: diversity,
        constraintValid,
        totalContribution: constraintValid ? total : 0,
      });
    }
    
    return contributions;
  }

  /**
   * Robustness: Persistence across resampling
   * Outputs that appear in high-density, stable modes contribute more
   */
  private calculateRobustnessContribution(
    embedding: number[],
    modes: DistributionMode[]
  ): number {
    if (modes.length === 0) return 0;
    
    // Find nearest mode
    let minDistance = Infinity;
    let nearestMode: DistributionMode | null = null;
    
    for (const mode of modes) {
      const distance = this.euclideanDistance(embedding, mode.center);
      if (distance < minDistance) {
        minDistance = distance;
        nearestMode = mode;
      }
    }
    
    if (!nearestMode) return 0;
    
    // Robustness = density × robustness × (1 - normalized distance)
    const normalizedDistance = Math.min(minDistance / (nearestMode.center.length || 1), 1);
    return nearestMode.density * nearestMode.robustness * (1 - normalizedDistance);
  }

  /**
   * Novelty: Exploration of low-density regions
   * Outputs far from existing modes contribute more
   */
  private calculateNoveltyContribution(
    embedding: number[],
    modes: DistributionMode[],
    allEmbeddings: Map<string, number[]>
  ): number {
    if (modes.length === 0) return 1; // First output is novel
    
    // Distance to nearest mode center
    const minDistance = Math.min(
      ...modes.map(m => this.euclideanDistance(embedding, m.center))
    );
    
    // Local density (k-nearest neighbors)
    const localDensity = this.estimateLocalDensity(embedding, allEmbeddings);
    
    // Novelty = distance × (1 - density)
    // Normalize distance to [0, 1]
    // For unit vectors, max Euclidean distance = 2 (opposite directions)
    const maxDistance = 2.0;
    const normalizedDistance = Math.min(minDistance / maxDistance, 1.0);
    
    return normalizedDistance * (1 - localDensity);
  }

  /**
   * Diversity: Expansion of semantic coverage
   * Outputs that increase average pairwise distance contribute more
   */
  private calculateDiversityContribution(
    embedding: number[],
    allEmbeddings: Map<string, number[]>
  ): number {
    if (allEmbeddings.size < 2) return 1;
    
    // Average distance to all other outputs
    const distances: number[] = [];
    for (const otherEmbedding of allEmbeddings.values()) {
      if (otherEmbedding !== embedding) {
        distances.push(this.euclideanDistance(embedding, otherEmbedding));
      }
    }
    
    if (distances.length === 0) return 0;
    
    const avgDistance = distances.reduce((a, b) => a + b, 0) / distances.length;
    
    // Normalize to [0, 1]
    // For unit vectors, max Euclidean distance = 2 (opposite directions)
    const maxDistance = 2.0;
    return Math.min(avgDistance / maxDistance, 1.0);
  }

  /**
   * Estimate local density using k-nearest neighbors
   */
  private estimateLocalDensity(
    embedding: number[],
    allEmbeddings: Map<string, number[]>,
    k: number = 5
  ): number {
    const distances: number[] = [];
    
    for (const otherEmbedding of allEmbeddings.values()) {
      if (otherEmbedding !== embedding) {
        distances.push(this.euclideanDistance(embedding, otherEmbedding));
      }
    }
    
    if (distances.length === 0) return 0;
    
    distances.sort((a, b) => a - b);
    const kNearest = distances.slice(0, Math.min(k, distances.length));
    const avgKNearest = kNearest.reduce((a, b) => a + b, 0) / kNearest.length;
    
    // Density = inverse of average distance to k-nearest
    // Normalize to [0, 1]
    return 1 / (1 + avgKNearest);
  }

  /**
   * Check constraints
   * 
   * Validates output against basic constraints:
   * - Non-empty output
   * - Valid data type
   * 
   * Extended validation can be added:
   * - JSON schema validation
   * - Safety/content filtering
   * - Format requirements
   * - Size limits
   */
  private async checkConstraints(output: MonteCarloOutput): Promise<boolean> {
    // Basic constraint checks
    if (!output.output) return false;
    if (typeof output.output === 'string' && output.output.length === 0) return false;
    if (typeof output.output === 'object' && Object.keys(output.output).length === 0) return false;
    
    // Additional checks can be added here:
    // - Schema validation against taskFormat.outputSchema
    // - Safety checks (content filtering, toxicity detection)
    // - Format validation (required fields, types)
    // - Size limits (max length, max tokens, etc.)
    
    return true;
  }

  /**
   * Sample outputs based on user preference
   * Safe → high-density, high-robustness
   * Novel → low-density but valid
   * Diverse → maximize semantic spread
   */
  async sampleByPreference(
    outputs: MonteCarloOutput[],
    contributions: Map<string, ContributionScore>,
    preference: UserPreference
  ): Promise<string[]> {
    const { alpha, beta, gamma } = preference.preferenceVector;
    
    // Calculate preference-weighted scores
    const weightedScores = outputs.map(output => {
      const contrib = contributions.get(output.outputId);
      if (!contrib || !contrib.constraintValid) {
        return { outputId: output.outputId, score: 0 };
      }
      
      const score = alpha * contrib.robustnessContribution +
                    beta * contrib.noveltyContribution +
                    gamma * contrib.diversityContribution;
      
      return { outputId: output.outputId, score };
    });
    
    // Filter and sort
    const validOutputs = weightedScores
      .filter(o => o.score > 0)
      .sort((a, b) => b.score - a.score);
    
    // Sample based on policy (implicit in preference vector)
    // If alpha is high → safe (robustness)
    // If beta is high → novel (novelty)
    // If gamma is high → diverse (diversity)
    
    // Return top-K outputs
    const k = Math.min(5, validOutputs.length);
    return validOutputs.slice(0, k).map(o => o.outputId);
  }
}