sublinear 0.2.0

const path = require('path');
const { ConvergenceDetector } = require('./convergence/convergence-detector');
const { MetricsReporter } = require('./convergence/metrics-reporter');
const { MatrixUtils } = require('./utils/matrix-utils');

// Import the existing WASM solver when available, fallback to stub
let WasmSolver;
let wasmAvailable = false;

try {
  // Try to import the WASM solver
  const wasmPath = path.join(__dirname, '../js/solver.js');
  // Note: This would need to be adjusted for actual WASM module loading in Node.js
  // For now, we'll create a compatible interface
  wasmAvailable = false; // Disable for CLI demo
} catch (error) {
  console.warn('WASM solver not available, using JavaScript fallback');
  wasmAvailable = false;
}

/**
 * JavaScript fallback solver for CLI and server use
 */
class JSSolver {
  constructor(config = {}) {
    this.config = {
      matrix: config.matrix,
      method: config.method || 'jacobi',
      tolerance: config.tolerance || 1e-10,
      maxIterations: config.maxIterations || 1000,
      enableVerification: config.enableVerification || false,
      verbose: config.verbose || false,
      ...config
    };

    this.initialized = false;
    this.currentSolution = null;
    this.currentIteration = 0;
    this.currentResidual = Infinity;
    this.converged = false;

    // Initialize convergence detection and metrics
    this.convergenceDetector = new ConvergenceDetector({
      tolerance: this.config.tolerance,
      maxIterations: this.config.maxIterations,
      relativeToleranceEnabled: this.config.relativeToleranceEnabled !== false
    });

    this.metricsReporter = new MetricsReporter({
      verbose: this.config.verbose,
      enableProfiling: true
    });
  }

  async initialize() {
    if (this.initialized) return;

    // Validate matrix and setup solver
    this.validateMatrix(this.config.matrix);
    this.initialized = true;
  }

  validateMatrix(matrix) {
    if (!matrix) {
      throw new Error('Matrix is required');
    }

    if (!matrix.rows || !matrix.cols) {
      throw new Error('Matrix must have rows and cols properties');
    }

    if (!matrix.data) {
      throw new Error('Matrix must have data property');
    }

    // Validate matrix format
    if (matrix.format === 'dense') {
      if (!Array.isArray(matrix.data) || matrix.data.length !== matrix.rows) {
        throw new Error('Dense matrix data must be array of rows');
      }
    } else if (matrix.format === 'coo') {
      if (!matrix.data.values || !matrix.data.rowIndices || !matrix.data.colIndices) {
        throw new Error('COO matrix must have values, rowIndices, and colIndices');
      }
    } else {
      throw new Error(`Unsupported matrix format: ${matrix.format}`);
    }

    // Enhanced diagonal validation
    const diagonalValidation = MatrixUtils.validateDiagonalElements(matrix);

    if (!diagonalValidation.valid) {
      const issues = [];

      if (diagonalValidation.missingDiagonals.length > 0) {
        issues.push(`Missing diagonal elements at rows: ${diagonalValidation.missingDiagonals.join(', ')}`);
      }

      if (diagonalValidation.smallDiagonals.length > 0) {
        const smallDiagInfo = diagonalValidation.smallDiagonals
          .map(d => `row ${d.index}: ${d.value}`)
          .join(', ');
        issues.push(`Near-zero diagonal elements: ${smallDiagInfo}`);
      }

      // Auto-fix option for diagonal issues
      if (this.config.autoFixMatrix !== false) {
        if (this.config.verbose) {
          console.warn('Matrix has diagonal issues, attempting auto-fix...');
          console.warn('Issues found:', issues.join('; '));
        }

        try {
          const fixResult = MatrixUtils.ensureDiagonalDominance(matrix, {
            strategy: this.config.diagonalStrategy || 'rowsum_plus_one',
            verbose: this.config.verbose
          });

          // Update the matrix in config with the fixed version
          this.config.matrix = fixResult.matrix;

          if (this.config.verbose) {
            console.log(`Applied ${fixResult.fixes.length} diagonal fixes`);
          }

          return; // Matrix is now valid
        } catch (fixError) {
          if (this.config.verbose) {
            console.error('Auto-fix failed:', fixError.message);
          }
        }
      }

      throw new Error(`Matrix validation failed: ${issues.join('; ')}. ` +
        `Set config.autoFixMatrix=true to enable automatic fixes.`);
    }

    // Check conditioning
    const conditioning = MatrixUtils.analyzeConditioning(matrix);

    if (this.config.verbose && conditioning.conditioningGrade !== 'A') {
      console.warn(`Matrix conditioning grade: ${conditioning.conditioningGrade}`);
      console.warn('Recommendations:', conditioning.recommendations.join(', '));

      if (!conditioning.isDiagonallyDominant) {
        console.warn(`Diagonal dominance ratio: ${conditioning.diagonalDominanceRatio.toFixed(2)} (should be ≤ 1.0)`);
      }
    }

    // Store conditioning info for solver optimization
    this.matrixConditioning = conditioning;

    // Check if matrix is symmetric (important for CG)
    this.isSymmetric = MatrixUtils.isSymmetric(matrix);

    if (this.config.verbose) {
      console.log(`Matrix symmetry: ${this.isSymmetric ? 'Yes' : 'No'}`);
      if (!this.isSymmetric && this.config.method === 'conjugate-gradient') {
        console.warn('Warning: Conjugate Gradient works best with symmetric positive definite matrices');
      }
    }
  }

  async solve(vector, options = {}) {
    await this.initialize();

    const onProgress = options.onProgress || (() => {});
    const signal = options.signal;

    // Reset convergence detector and start metrics collection
    this.convergenceDetector.reset();
    this.convergenceDetector.initialize(vector);

    this.metricsReporter.startSolve(this.config, {
      rows: this.config.matrix.rows,
      cols: this.config.matrix.cols,
      format: this.config.matrix.format || 'dense'
    });

    // Initialize solution vector
    this.currentSolution = new Array(vector.length).fill(0);
    this.currentIteration = 0;
    this.currentResidual = Infinity;
    this.converged = false;

    const startTime = Date.now();

    try {
      // Choose solver method
      switch (this.config.method) {
        case 'jacobi':
          await this.solveJacobi(vector, onProgress, signal);
          break;
        case 'gauss-seidel':
        case 'gauss_seidel':
          await this.solveGaussSeidel(vector, onProgress, signal);
          break;
        case 'conjugate-gradient':
        case 'cg':
          await this.solveConjugateGradient(vector, onProgress, signal);
          break;
        case 'hybrid':
        case 'adaptive':
        default:
          await this.solveAdaptive(vector, onProgress, signal);
          break;
      }

      // Generate final report
      const report = this.metricsReporter.finalizeSolve(this.convergenceDetector, this.currentSolution);
      const elapsed = Date.now() - startTime;

      return {
        values: this.currentSolution,
        iterations: this.currentIteration,
        residual: this.currentResidual,
        converged: this.converged,
        solveTime: elapsed,
        method: this.config.method,
        memoryUsage: this.getMemoryUsage(),

        // Enhanced metrics
        convergenceReport: report,
        convergenceRate: report.convergence.convergenceRatePercent,
        reductionFactor: report.convergence.reductionFactor,
        performanceGrade: report.performance.grade
      };

    } catch (error) {
      throw new Error(`Solve failed: ${error.message}`);
    }
  }

  async *streamSolve(vector) {
    await this.initialize();

    // Initialize solution vector
    this.currentSolution = new Array(vector.length).fill(0);
    this.currentIteration = 0;
    this.currentResidual = Infinity;
    this.converged = false;

    const startTime = Date.now();

    try {
      // Stream solver based on method
      switch (this.config.method) {
        case 'jacobi':
          yield* this.streamJacobi(vector);
          break;
        case 'gauss-seidel':
        case 'gauss_seidel':
          yield* this.streamGaussSeidel(vector);
          break;
        case 'conjugate-gradient':
        case 'cg':
          yield* this.streamConjugateGradient(vector);
          break;
        case 'hybrid':
        case 'adaptive':
        default:
          yield* this.streamAdaptive(vector);
          break;
      }

    } catch (error) {
      yield {
        error: error.message,
        iteration: this.currentIteration,
        timestamp: new Date().toISOString()
      };
    }
  }

  async *streamJacobi(vector) {
    const matrix = this.config.matrix;
    const n = vector.length;
    let x = new Array(n).fill(0);
    let xNew = new Array(n).fill(0);

    for (let iter = 0; iter < this.config.maxIterations; iter++) {
      this.currentIteration = iter;

      // Jacobi iteration
      for (let i = 0; i < n; i++) {
        let sum = 0;
        let diagonal = 0;

        if (matrix.format === 'dense') {
          for (let j = 0; j < n; j++) {
            if (i !== j) {
              sum += matrix.data[i][j] * x[j];
            } else {
              diagonal = matrix.data[i][j];
            }
          }
        } else if (matrix.format === 'coo') {
          // Coordinate format
          for (let k = 0; k < matrix.data.values.length; k++) {
            const row = matrix.data.rowIndices[k];
            const col = matrix.data.colIndices[k];
            const val = matrix.data.values[k];

            if (row === i) {
              if (col !== i) {
                sum += val * x[col];
              } else {
                diagonal = val;
              }
            }
          }
        }

        if (Math.abs(diagonal) < 1e-14) {
          throw new Error(`Zero diagonal element at position ${i}`);
        }

        xNew[i] = (vector[i] - sum) / diagonal;
      }

      // Update solution
      x = [...xNew];
      this.currentSolution = x;

      // Use convergence detector for proper metrics
      const convergenceMetrics = this.convergenceDetector.update(matrix, x, vector);
      this.currentIteration = convergenceMetrics.iteration;
      this.currentResidual = convergenceMetrics.relativeResidualNorm;
      this.converged = convergenceMetrics.isConverged;

      // Record metrics
      const iterationMetrics = this.metricsReporter.recordIteration(convergenceMetrics);

      // Yield progress with enhanced metrics
      yield {
        iteration: iter,
        residual: convergenceMetrics.relativeResidualNorm,
        residualNorm: convergenceMetrics.residualNorm,
        convergenceRate: convergenceMetrics.convergenceRate,
        convergenceRatePercent: (1 - convergenceMetrics.convergenceRate) * 100,
        reductionFactor: convergenceMetrics.reductionFactor,
        memoryUsage: this.getMemoryUsage(),
        converged: this.converged,
        shouldStop: convergenceMetrics.shouldStop,
        estimatedIterationsRemaining: convergenceMetrics.estimatedIterationsRemaining,
        solution: this.converged ? x : undefined,
        verified: this.config.enableVerification ? await this.verify(x, vector) : undefined
      };

      // Early stopping when convergence criteria met
      if (convergenceMetrics.shouldStop) {
        break;
      }

      // Allow other operations to run
      await new Promise(resolve => setImmediate(resolve));
    }
  }

  async *streamGaussSeidel(vector) {
    const matrix = this.config.matrix;
    const n = vector.length;
    let x = new Array(n).fill(0);

    for (let iter = 0; iter < this.config.maxIterations; iter++) {
      this.currentIteration = iter;

      // Gauss-Seidel iteration
      for (let i = 0; i < n; i++) {
        let sum = 0;
        let diagonal = 0;

        if (matrix.format === 'dense') {
          for (let j = 0; j < n; j++) {
            if (i !== j) {
              sum += matrix.data[i][j] * x[j];
            } else {
              diagonal = matrix.data[i][j];
            }
          }
        } else if (matrix.format === 'coo') {
          for (let k = 0; k < matrix.data.values.length; k++) {
            const row = matrix.data.rowIndices[k];
            const col = matrix.data.colIndices[k];
            const val = matrix.data.values[k];

            if (row === i) {
              if (col !== i) {
                sum += val * x[col];
              } else {
                diagonal = val;
              }
            }
          }
        }

        if (Math.abs(diagonal) < 1e-14) {
          throw new Error(`Zero diagonal element at position ${i}`);
        }

        x[i] = (vector[i] - sum) / diagonal;
      }

      // Compute residual
      const residual = this.computeResidual(matrix, x, vector);
      this.currentResidual = this.vectorNorm(residual);

      // Convergence check
      this.converged = this.currentResidual < this.config.tolerance;

      // Update solution
      this.currentSolution = [...x];

      // Yield progress
      yield {
        iteration: iter,
        residual: this.currentResidual,
        convergenceRate: iter > 0 ? this.currentResidual / this.previousResidual : 1.0,
        memoryUsage: this.getMemoryUsage(),
        converged: this.converged,
        solution: this.converged ? x : undefined,
        verified: this.config.enableVerification ? await this.verify(x, vector) : undefined
      };

      this.previousResidual = this.currentResidual;

      if (this.converged) {
        break;
      }

      await new Promise(resolve => setImmediate(resolve));
    }
  }

  async *streamConjugateGradient(vector) {
    const matrix = this.config.matrix;
    const n = vector.length;
    let x = new Array(n).fill(0);

    // CG initialization: r = b - Ax (initial residual)
    let r = [...vector]; // Since x starts at 0, r = b - A*0 = b
    let p = [...r];      // Initial search direction
    let rsold = this.dotProduct(r, r);

    for (let iter = 0; iter < this.config.maxIterations; iter++) {
      this.currentIteration = iter;

      // Compute Ap (matrix-vector product)
      const Ap = this.multiplyMatrixVector(matrix, p);
      const pAp = this.dotProduct(p, Ap);

      // Check for non-positive curvature (matrix not positive definite)
      if (pAp <= 1e-16) {
        if (this.config.verbose) {
          console.warn(`CG: Non-positive curvature detected at iteration ${iter}, switching to steepest descent`);
        }
        // Fall back to steepest descent step
        const rAr = this.dotProduct(r, this.multiplyMatrixVector(matrix, r));
        if (rAr > 1e-16) {
          const alpha = rsold / rAr;
          for (let i = 0; i < n; i++) {
            x[i] += alpha * r[i];
          }
        }
        break;
      }

      // CG step size
      const alpha = rsold / pAp;

      // Update solution: x = x + alpha * p
      for (let i = 0; i < n; i++) {
        x[i] += alpha * p[i];
      }

      // Update residual: r = r - alpha * Ap
      for (let i = 0; i < n; i++) {
        r[i] -= alpha * Ap[i];
      }

      const rsnew = this.dotProduct(r, r);
      this.currentResidual = Math.sqrt(rsnew);

      // Update solution
      this.currentSolution = [...x];

      // Use convergence detector for proper metrics
      const convergenceMetrics = this.convergenceDetector.update(matrix, x, vector);
      this.currentIteration = convergenceMetrics.iteration;
      this.currentResidual = convergenceMetrics.relativeResidualNorm;
      this.converged = convergenceMetrics.isConverged;

      // Record metrics
      const iterationMetrics = this.metricsReporter.recordIteration(convergenceMetrics);

      // Yield progress with enhanced metrics
      yield {
        iteration: iter,
        residual: convergenceMetrics.relativeResidualNorm,
        residualNorm: convergenceMetrics.residualNorm,
        convergenceRate: convergenceMetrics.convergenceRate,
        convergenceRatePercent: (1 - convergenceMetrics.convergenceRate) * 100,
        reductionFactor: convergenceMetrics.reductionFactor,
        memoryUsage: this.getMemoryUsage(),
        converged: this.converged,
        shouldStop: convergenceMetrics.shouldStop,
        estimatedIterationsRemaining: convergenceMetrics.estimatedIterationsRemaining,
        solution: this.converged ? x : undefined,
        verified: this.config.enableVerification ? await this.verify(x, vector) : undefined
      };

      // Early stopping when convergence criteria met
      if (convergenceMetrics.shouldStop) {
        break;
      }

      // Check for stagnation
      if (rsnew > rsold * 0.999) {
        if (this.config.verbose) {
          console.warn(`CG: Convergence stagnation detected at iteration ${iter}`);
        }
      }

      // Update search direction: p = r + beta * p
      const beta = rsnew / rsold;
      for (let i = 0; i < n; i++) {
        p[i] = r[i] + beta * p[i];
      }

      rsold = rsnew;
      await new Promise(resolve => setImmediate(resolve));
    }
  }

  async *streamAdaptive(vector) {
    // Start with Jacobi, switch to CG if convergence is slow
    let method = 'jacobi';
    let slowConvergence = false;
    let previousResiduals = [];

    const jacobiSolver = new JSSolver({
      ...this.config,
      method: 'jacobi'
    });

    await jacobiSolver.initialize();

    for await (const update of jacobiSolver.streamSolve(vector)) {
      previousResiduals.push(update.residual);

      // Check for slow convergence after 50 iterations
      if (update.iteration > 50 && !slowConvergence) {
        const recent = previousResiduals.slice(-10);
        const improvement = recent[0] / recent[recent.length - 1];

        if (improvement < 1.1) { // Less than 10% improvement
          slowConvergence = true;
          method = 'conjugate-gradient';

          // Switch to CG
          const cgSolver = new JSSolver({
            ...this.config,
            method: 'conjugate-gradient'
          });

          await cgSolver.initialize();

          // Continue with CG from current solution
          const remainingVector = vector; // In practice, would adjust for current progress

          for await (const cgUpdate of cgSolver.streamSolve(remainingVector)) {
            yield {
              ...cgUpdate,
              iteration: update.iteration + cgUpdate.iteration,
              method: 'adaptive-cg'
            };

            if (cgUpdate.converged) {
              return;
            }
          }
          return;
        }
      }

      yield {
        ...update,
        method: 'adaptive-jacobi'
      };

      if (update.converged) {
        return;
      }
    }
  }

  // Non-streaming versions
  async solveJacobi(vector, onProgress, signal) {
    for await (const update of this.streamJacobi(vector)) {
      if (signal && signal.aborted) {
        throw new Error('Solve aborted');
      }

      onProgress(update);

      if (update.converged) {
        break;
      }
    }
  }

  async solveGaussSeidel(vector, onProgress, signal) {
    for await (const update of this.streamGaussSeidel(vector)) {
      if (signal && signal.aborted) {
        throw new Error('Solve aborted');
      }

      onProgress(update);

      if (update.converged) {
        break;
      }
    }
  }

  async solveConjugateGradient(vector, onProgress, signal) {
    for await (const update of this.streamConjugateGradient(vector)) {
      if (signal && signal.aborted) {
        throw new Error('Solve aborted');
      }

      onProgress(update);

      if (update.converged) {
        break;
      }
    }
  }

  async solveAdaptive(vector, onProgress, signal) {
    for await (const update of this.streamAdaptive(vector)) {
      if (signal && signal.aborted) {
        throw new Error('Solve aborted');
      }

      onProgress(update);

      if (update.converged) {
        break;
      }
    }
  }

  // Utility methods
  computeResidual(matrix, x, b) {
    const Ax = this.multiplyMatrixVector(matrix, x);
    return Ax.map((val, i) => b[i] - val);
  }

  multiplyMatrixVector(matrix, vector) {
    const result = new Array(matrix.rows).fill(0);

    if (matrix.format === 'dense') {
      for (let i = 0; i < matrix.rows; i++) {
        for (let j = 0; j < matrix.cols; j++) {
          result[i] += matrix.data[i][j] * vector[j];
        }
      }
    } else if (matrix.format === 'coo') {
      for (let k = 0; k < matrix.data.values.length; k++) {
        const row = matrix.data.rowIndices[k];
        const col = matrix.data.colIndices[k];
        const val = matrix.data.values[k];
        result[row] += val * vector[col];
      }
    }

    return result;
  }

  vectorNorm(vector) {
    return Math.sqrt(vector.reduce((sum, val) => sum + val * val, 0));
  }

  dotProduct(a, b) {
    return a.reduce((sum, val, i) => sum + val * b[i], 0);
  }

  async verify(solution, vector) {
    if (!this.config.enableVerification) {
      return { verified: true };
    }

    // Simple verification: check residual
    const residual = this.computeResidual(this.config.matrix, solution, vector);
    const residualNorm = this.vectorNorm(residual);

    return {
      verified: residualNorm < this.config.tolerance * 10,
      residualNorm,
      maxError: Math.max(...residual.map(Math.abs))
    };
  }

  getMemoryUsage() {
    const memUsage = process.memoryUsage();
    return Math.round(memUsage.heapUsed / 1024 / 1024); // MB
  }

  stop() {
    // For stopping streaming solves
    this.stopped = true;
  }
}

/**
 * Factory function to create solver instances
 */
async function createSolver(config = {}) {
  const solver = wasmAvailable
    ? new WasmSolver(config)
    : new JSSolver(config);

  await solver.initialize();
  return solver;
}

module.exports = {
  createSolver,
  JSSolver,
  WasmSolver: wasmAvailable ? WasmSolver : JSSolver
};