// =============================================================
// СТАТИСТИЧЕСКОЕ СОГЛАШЕНИЕ (LinearAlgebraML)
//
// Используется дисперсия генеральной совокупности (деление на n),
// как принято в численных методах и алгоритмах машинного обучения.
//
// Исключение: Matrix.PCA использует (n-1) для ковариационной матрицы
// согласно статистической политике MLABC (пункт 3).
//
// См. статистическую политику в модуле MLABC.
// =============================================================

/// Линейная алгебра для алгоритмов машинного обучения.
///
/// Содержит типы Vector и Matrix и численные методы,
/// используемые в моделях ML:
/// решение систем линейных уравнений, least squares,
/// ridge-регрессия и специализированные методы
/// для симметричных положительно определённых матриц.
///
/// Модуль предназначен для внутреннего использования
/// в ML-алгоритмах и оптимизирован для задач обучения моделей.
unit LinearAlgebraML;

interface

type
  Vector = class
  private
    fdata: array of real;
    static procedure CheckSameLength(const a, b: Vector);
    static procedure CheckNonEmpty(const v: Vector); 
    procedure SetData(i: integer; value: real) := fdata[i] := value;
  public
    /// Массив элементов вектора.
    property Data: array of real read fdata;
    /// Длина вектора.
    property Length: integer read fdata.Length;
    /// Доступ к элементу вектора по индексу.
    property Item[i: integer]: real read fdata[i] write SetData; default;
    
    /// Создаёт вектор заданной длины, заполненный нулями.
    constructor Create(n: integer);
    /// Создаёт вектор из массива вещественных значений.
    constructor Create(values: array of real);
    /// Создаёт вектор из массива целых значений.
    constructor Create(values: array of integer);
    
    /// Возвращает копию вектора.
    function Clone: Vector;
    
    /// Нормирует текущий вектор по евклидовой норме и возвращает его.
    function Normalize: Vector;
    /// Возвращает нормированную копию вектора.
    function Normalized: Vector;
    
    /// Возвращает копию данных в виде массива.
    function ToArray: array of real;
    /// Преобразует элементы вектора в массив целых чисел.
    function ToIntArray: array of integer;
    
    static function operator +(a, b: Vector): Vector;
    static function operator -(a, b: Vector): Vector;
    
    static function operator *(alpha: real; v: Vector): Vector;
    static function operator *(v: Vector; alpha: real): Vector;
    static function operator /(v: Vector; alpha: real): Vector;
    
    static function operator +=(a: Vector; b: Vector): Vector;
    static function operator -=(a: Vector; b: Vector): Vector;
    static function operator *=(a: Vector; alpha: real): Vector;
    
    static function operator +(v: Vector; c: real): Vector;
    static function operator +(c: real; v: Vector): Vector;
    
    static function operator implicit(a: array of real): Vector := new Vector(a);
    static function operator implicit(a: array of integer): Vector := new Vector(a);

    // ---------- Векторные функции ----------
    /// Применить функцию ко всем элементам вектора
    function Apply(f: real -> real): Vector;
    /// Возвращает вектор из квадратных корней элементов.
    function Sqrt: Vector;
    /// Возвращает вектор из экспонент элементов.
    function Exp: Vector := Apply(PABCSystem.Exp);
    /// Возвращает вектор из натуральных логарифмов элементов.
    function Ln: Vector := Apply(PABCSystem.Ln);
    /// Возвращает вектор из модулей элементов.
    function Abs: Vector := Apply(PABCSystem.Abs);    
    
    // ---------- Основные методы ----------
    /// Возвращает сумму элементов вектора.
    function Sum: real;
    /// Среднее. Синоним Mean
    function Average: real;
    /// Среднее
    function Mean: real;
    /// Возвращает квадрат евклидовой нормы вектора.
    function Norm2: real;
    /// Возвращает евклидову норму вектора.
    function Norm: real;
    /// Возвращает максимальный элемент вектора.
    function Max: real;
    /// Возвращает минимальный элемент вектора.
    function Min: real;
    /// Скалярное произведение
    function Dot(b: Vector): real;
    
    // ---------- Сервисные методы ----------
    /// Возвращает строковое представление вектора.
    function ToString: string; override := $'{fdata.Select(x -> x.ToString(''G3''))}';
    /// Возвращает строковое представление вектора с заданной точностью.
    function ToString(digits: integer): string := $'{fdata.Select(x -> x.ToString(''G''+digits))}';
    /// Печатает вектор без перевода строки.
    procedure Print := fdata.Print;
    /// Печатает вектор с переводом строки.
    procedure Println := fdata.Println;
    
    /// Возвращает подмножество элементов по заданным индексам.
    function SubvectorBy(indices: array of integer): Vector;
  end;
  
  Matrix = class
  private
    fdata: array[,] of real;

    static procedure CheckSameSize(A, B: Matrix);
    static procedure CheckMulSize(A, B: Matrix);
    static procedure CheckVecSize(A: Matrix; x: Vector);
    
    procedure SetData(i, j: integer; value: real) := fdata[i, j] := value;
  public
    /// Число строк матрицы.
    property RowCount: integer read fdata.RowCount;
    /// Число столбцов матрицы.
    property ColCount: integer read fdata.ColCount;
    /// Данные матрицы в виде двумерного массива.
    property Data: array[,] of real read fdata;
    
    /// Доступ к элементу матрицы по индексам строки и столбца.
    property Item[i, j: integer]: real
      read fdata[i, j] write SetData; default;
    
    /// Создаёт матрицу заданного размера, заполненную нулями.
    constructor Create(r, c: integer);
    /// Создаёт матрицу из двумерного массива вещественных значений.
    constructor Create(values: array[,] of real);
    
    constructor Create(values: array of Vector);
    
    /// Возвращает копию данных в виде двумерного массива.
    function ToArray2D: array[,] of real;
    /// Возвращает строку матрицы в виде массива.
    function RowToArray(r: integer): array of real;
    
    /// Возвращает копию матрицы.
    function Clone: Matrix;
    
    /// Возвращает строку матрицы в виде массива.
    function Row(i: integer): array of real := fdata.Row(i);
    /// Возвращает столбец матрицы в виде массива.
    function Col(j: integer): array of real := fdata.Col(j);

    /// Возвращает суммы по всем столбцам.
    function ColumnSums: Vector;
    /// Возвращает суммы по всем строкам.
    function RowSums: Vector;
    /// Возвращает средние значения по всем столбцам.
    function ColumnMeans: Vector;
    /// Возвращает средние значения по всем строкам.
    function RowMeans: Vector;
    /// Возвращает дисперсии по всем столбцам.
    function ColumnVariances: Vector;
    /// Возвращает дисперсии по всем строкам.
    function RowVariances: Vector;
    /// Возвращает стандартные отклонения по всем столбцам.
    function ColumnStd: Vector;
    /// Возвращает стандартные отклонения по всем строкам.
    function RowStd: Vector;
    /// Возвращает минимумы по всем столбцам.
    function ColumnMins: Vector;
    /// Возвращает максимумы по всем столбцам.
    function ColumnMaxs: Vector;
    /// Возвращает минимумы по всем строкам.
    function RowMins: Vector;
    /// Возвращает максимумы по всем строкам.
    function RowMaxs: Vector;
    
    /// Возвращает сумму элементов строки.
    function RowSum(i: integer): real;
    /// Возвращает среднее значение элементов строки.
    function RowMean(i: integer): real;
    /// Возвращает дисперсию элементов строки.
    function RowVariance(i: integer): real;
    /// Возвращает стандартное отклонение элементов строки.
    function RowStd(i: integer): real;
    /// Возвращает минимальный элемент строки.
    function RowMin(i: integer): real;
    /// Возвращает максимальный элемент строки.
    function RowMax(i: integer): real;
    /// Возвращает индекс минимального элемента строки.
    function RowArgMin(i: integer): integer;
    /// Возвращает индекс максимального элемента строки.
    function RowArgMax(i: integer): integer;
    
    /// Возвращает сумму элементов столбца.
    function ColumnSum(j: integer): real;
    /// Возвращает среднее значение элементов столбца.
    function ColumnMean(j: integer): real;
    /// Возвращает дисперсию элементов столбца.
    function ColumnVariance(j: integer): real;
    /// Возвращает стандартное отклонение элементов столбца.
    function ColumnStd(j: integer): real;
    /// Возвращает минимальный элемент столбца.
    function ColumnMin(j: integer): real;
    /// Возвращает максимальный элемент столбца.
    function ColumnMax(j: integer): real;
    /// Возвращает индекс минимального элемента столбца.
    function ColumnArgMin(j: integer): integer;
    /// Возвращает индекс максимального элемента столбца.
    function ColumnArgMax(j: integer): integer;
    
    /// Возвращает норму Фробениуса матрицы.
    function FrobeniusNorm: real;
    /// Прибавляет lambda к диагональным элементам матрицы.
    procedure AddScaledIdentity(lambda: real);

    
    static function operator implicit(a: array [,] of real): Matrix := new Matrix(a);
    static function operator implicit(a: array of Vector): Matrix := new Matrix(a);
    static function operator implicit(a: array of array of real): Matrix := new Matrix(Matr(a));
    static function operator implicit(a: array of array of integer): Matrix := new Matrix(Matr(a.ConvertAll(x -> x.ConvertAll(y -> real(y)))));
    
    // ---------- value operators ----------
    static function operator +(A, B: Matrix): Matrix;
    static function operator -(A, B: Matrix): Matrix;
    
    static function operator *(A: Matrix; x: Vector): Vector;
    static function operator *(A, B: Matrix): Matrix;
    
    // ---------- in-place operators ----------
    static function operator +=(A, B: Matrix): Matrix;
    static function operator -=(A, B: Matrix): Matrix;
    static function operator *=(A: Matrix; alpha: real): Matrix;
    
    // ---------- Основные методы ----------
    
    /// Возвращает транспонированную матрицу Aᵀ
    function Transpose: Matrix;
    /// Проверяет, является ли матрица симметричной с заданным допуском.
    function IsSymmetric(tol: real := 1e-12): boolean;
    /// Вычисляет собственные значения и собственные векторы вещественной квадратной симметричной матрицы.
    ///
    /// Возвращает:
    /// - values  — вектор длины n, содержащий собственные значения,
    ///             отсортированные по убыванию;
    /// - vectors — матрицу n × n, столбцы которой являются
    ///             ортонормированными собственными векторами.
    ///
    /// Выполняется разложение:
    ///     A = V * diag(values) * Vᵀ
    ///
    /// Метод основан на итерациях Якоби.    
    function EigenSymmetric(tol: real := 1e-12; maxIter: integer := 100): (Vector, Matrix);
    /// Выполняет анализ главных компонент (PCA) по матрице данных.
    ///
    /// Строки интерпретируются как объекты, столбцы — как признаки.
    /// Параметр k задаёт число главных компонент (1 ≤ k ≤ Cols).
    ///
    /// Возвращает:
    /// - components — матрицу Cols × k главных компонент;
    /// - variances  — соответствующие дисперсии.
    ///
    /// Компоненты ортонормированы и отсортированы по убыванию дисперсии.
    function PCA(k: integer): (Matrix, Vector);
    
    // ---------- Сервисные методы ----------
    /// Возвращает строковое представление матрицы.
    function ToString: string; override := $'{fdata}';
    /// Печатает матрицу без перевода строки.
    procedure Print := fdata.Print;
    /// Печатает матрицу с переводом строки.
    procedure Println := fdata.Println;
    
    /// Возвращает строку матрицы в виде вектора.
    function GetRow(i: integer): Vector;
    /// Возвращает столбец матрицы в виде вектора.
    function GetCol(j: integer): Vector;
    
    /// Возвращает матрицу, составленную из строк с заданными индексами.
    function TakeRows(indices: array of integer): Matrix;

    // ---------- Статические методы ----------
    /// Возвращает единичную матрицу размера n
    static function Identity(n: integer): Matrix;
    /// Возвращает внешнее произведение двух векторов
    static function OuterProduct(a, b: Vector): Matrix;
    
  end;

/// Решает систему линейных уравнений A * x = b с помощью LU-разложения с частичным выбором главного элемента.
/// A должна быть квадратной матрицей. 
///
/// Временная сложность: O(n³).
///
/// Вызывает исключение, если A вырождена или размеры не согласованы
function Solve(A: Matrix; b: Vector): Vector;

/// Решает систему линейных уравнений A * x = b, предполагая, что матрица A является
/// симметричной положительно определённой (SPD).
/// Используется разложение Холецкого (A = L * L^T).
///
/// Временная сложность: O(n³ / 3).
///
/// Вызывает исключение, если матрица A не симметричная положительно определенная
/// или если размеры A и b не согласованы.
function SolveSPD(A: Matrix; b: Vector): Vector;

/// Решает систему линейных уравнений A * x = b.
///
/// Метод автоматически выбирается для обеспечения
/// максимальной производительности и численной устойчивости.
///
/// Вызывает исключение, если система неразрешима
/// или размеры A и b не согласованы.
function SolveAuto(A: Matrix; b: Vector): Vector;


/// Решает систему A * x ≈ b методом ridge-регрессии с подавлением неустойчивых направлений.
/// 
/// Матрица A имеет размер m × n: m — число уравнений (строк A), n — число неизвестных (столбцов A).
///
/// Метод используется, когда обычное решение (lambda = 0) даёт слишком большие коэффициенты из-за  
/// сильной зависимости между столбцами матрицы A или из-за шума в данных.
///
/// Параметр lambda задаёт степень регуляризации:
/// - lambda = 0
///     эквивалентно обычному решению по методу наименьших квадратов;
/// - lambda ≈ 1e-6 .. 1e-3
///     слабая регуляризация;
/// - lambda ≈ 1e-2 .. 1
///     умеренная регуляризация (типично для практических и ML-задач);
/// - lambda > 1
///     сильная регуляризация, коэффициенты заметно уменьшаются.
///
/// Если данные хорошо обусловлены, используйте lambda = 0
/// При шумных или коррелированных данных часто подходят значения 0.01 .. 0.1
/// Увеличивайте lambda, если коэффициенты x становятся слишком большими.
///
/// Требуется, чтобы число уравнений было не меньше числа неизвестных (m ≥ n).
/// Входные A и b изменяются. Возвращаемый вектор x имеет длину n.
function SolveRidge(A: Matrix; b: Vector; lambda: real := 0): Vector;

/// Решает задачу наименьших квадратов для системы A * x ≈ b.
///
/// Используется QR-разложение матрицы A.
/// Подходит для переопределённых систем (число строк больше числа столбцов).
///
/// Временная сложность: O(m * n²), где
/// m — число строк A,
/// n — число столбцов A.
///
/// Вызывает исключение, если размеры A и b не согласованы.
function SolveLeastSquaresQR(A: Matrix; b: Vector): Vector;

implementation

uses MLExceptions;

const
  ER_VECTOR_LENGTH_NEGATIVE =
    'Длина вектора должна быть неотрицательной!!Vector length must be non-negative';
  ER_VALUES_NULL =
    'values не может быть nil!!values cannot be nil';
  ER_VALUES_EMPTY =
    'values не может быть пустым!!values cannot be empty';
  ER_VECTOR_LENGTH_MISMATCH =
    'Несоответствие длины векторов: {0} и {1}!!Vector length mismatch: {0} vs {1}';
  ER_VECTOR_NULL =
    'Вектор не может быть nil!!vector cannot be nil';
  ER_MATRIX_ROWS_HAVE_DIFFERENT_LENGTHS =
    'Строки матрицы имеют разную длину!!matrix rows have different lengths';  
  ER_VECTOR_EMPTY =
    'Вектор пуст!!Vector is empty';
  ER_VECTOR_DIVIDE_BY_ZERO =
    'Деление на ноль при делении вектора на скаляр!!Division by zero in Vector / scalar';
  ER_VECTOR_TO_INT_NON_INTEGER =
    'Вектор содержит нецелое значение на позиции {0}: {1}!!Vector contains non-integer value at index {0}: {1}';
  ER_MATRIX_SIZE_NEGATIVE =
    'Размеры матрицы должны быть неотрицательными!!Matrix size must be non-negative';
  ER_MATRIX_SIZE_MISMATCH =
    'Несоответствие размеров матриц: {0}x{1} и {2}x{3}!!' +
    'Matrix size mismatch: {0}x{1} vs {2}x{3}';
  ER_MATRIX_MUL_SIZE_MISMATCH =
    'Несоответствие размеров при умножении матриц: {0}x{1} * {2}x{3}!!' +
    'Matrix multiply size mismatch: {0}x{1} * {2}x{3}';
  ER_MATRIX_VECTOR_SIZE_MISMATCH =
    'Несоответствие размеров при умножении матрицы на вектор: {0}x{1} * {2}!!' +
    'Matrix-vector size mismatch: {0}x{1} * {2}';
  ER_ROW_INDEX_OUT_OF_RANGE =
    'Индекс строки {0} вне диапазона [0..{1})!!Row index {0} out of range [0..{1})';
  ER_COL_INDEX_OUT_OF_RANGE =
    'Индекс столбца {0} вне диапазона [0..{1})!!Column index {0} out of range [0..{1})';
  ER_MATRIX_NOT_SQUARE =
    'Матрица должна быть квадратной!!Matrix must be square';
  ER_MATRIX_NOT_SYMMETRIC =
    'Матрица должна быть симметричной!!Matrix must be symmetric';
  ER_PCA_K_INVALID =
    'k должно быть >= 1!!k must be >= 1';
  ER_PCA_K_TOO_LARGE =
    'k не может превышать число признаков!!k cannot exceed number of features';
  ER_PCA_NEED_TWO_SAMPLES =
    'Требуется как минимум два объекта!!At least two samples required';
  ER_CHOLESKY_NOT_SQUARE =
    'Матрица должна быть квадратной для Cholesky!!Matrix must be square for Cholesky';
  ER_MATRIX_NOT_SPD =
    'Матрица не является положительно определённой (SPD)!!Matrix is not SPD';
  ER_MATRIX_SINGULAR =
    'Матрица вырождена!!Matrix is singular';
  ER_VECTOR_SIZE_MISMATCH =
    'Несоответствие длины вектора: {0} и {1}!!Vector size mismatch: {0} and {1}';
  ER_QR_REQUIRES_M_GE_N =
    'Для QR-разложения требуется m >= n!!QR decomposition requires m >= n';
  ER_SINGULAR_MATRIX =
    'Матрица вырождена или плохо обусловлена!!Matrix is singular or ill-conditioned';
  ER_EMPTY_MATRIX =
    'Матрица пуста!!Matrix is empty';    
  ER_EIGEN_NOT_CONVERGED =
    'EigenSymmetric: не сошлось за {0} итераций (off={1}, tol={2})!!EigenSymmetric did not converge in {0} iterations (off={1}, tol={2})';    
  
type
  MLNotSPDException = class(MLException);
    
//-----------------------------
//           Vector
//-----------------------------

constructor Vector.Create(n: integer);
begin
  if n < 0 then
    ArgumentOutOfRangeError(ER_VECTOR_LENGTH_NEGATIVE);
  fdata := new real[n];
end;

constructor Vector.Create(values: array of real);
begin
  if values = nil then
    ArgumentNullError(ER_VALUES_NULL);
  fdata := Copy(values);
end;

constructor Vector.Create(values: array of integer);
begin
  if values = nil then
    ArgumentNullError(ER_VALUES_NULL);
  fdata := values.Select(x -> real(x)).ToArray;
end;

function Vector.ToArray: array of real;
begin
  Result := Copy(fdata);
end;

function Vector.ToIntArray: array of integer;
begin
  Result := new integer[Length];
  
  for var i := 0 to Length - 1 do
  begin
    var x := Data[i];
    var ix := Round(x);
    
    if PABCSystem.Abs(x - ix) > 1e-12 then
      ArgumentError(ER_VECTOR_TO_INT_NON_INTEGER, i, x);
    
    Result[i] := ix;
  end;
end;

function Vector.Clone: Vector;
begin
  Result := new Vector(fdata);
end;

function Vector.Normalize: Vector;
begin
  var sum := 0.0;
  for var i := 0 to Length - 1 do
    sum += Self[i];

  if sum > 0 then
    for var i := 0 to Length - 1 do
      Self[i] /= sum;

  Result := Self;
end;

function Vector.Normalized: Vector;
begin
  Result := Self.Clone;
  Result.Normalize;
end;

static procedure Vector.CheckSameLength(a, b: Vector);
begin
  if a.Length <> b.Length then
    ArgumentError(ER_VECTOR_LENGTH_MISMATCH, a.Length, b.Length);
end;

static procedure Vector.CheckNonEmpty(v: Vector);
begin
  if v.Length = 0 then
    ArgumentError(ER_VECTOR_EMPTY);
end;

static function Vector.operator +(a, b: Vector): Vector;
begin
  CheckSameLength(a, b);
  Result := new Vector(a.Length);
  for var i := 0 to a.Length - 1 do
    Result.fdata[i] := a.fdata[i] + b.fdata[i];
end;

static function Vector.operator -(a, b: Vector): Vector;
begin
  CheckSameLength(a, b);
  Result := new Vector(a.Length);
  for var i := 0 to a.Length - 1 do
    Result.fdata[i] := a.fdata[i] - b.fdata[i];
end;

static function Vector.operator *(v: Vector; alpha: real): Vector;
begin
  Result := new Vector(v.Length);
  for var i := 0 to v.Length - 1 do
    Result.fdata[i] := alpha * v.fdata[i];
end;

static function Vector.operator *(alpha: real; v: Vector): Vector;
begin
  Result := v * alpha;
end;

static function Vector.operator /(v: Vector; alpha: real): Vector;
begin
  if alpha = 0.0 then
    raise new System.DivideByZeroException(GetTranslation(ER_VECTOR_DIVIDE_BY_ZERO));
  
  Result := new Vector(v.Length);
  var inv := 1.0 / alpha;
  for var i := 0 to v.Length - 1 do
    Result.fdata[i] := v.fdata[i] * inv;
end;

static function Vector.operator +=(a, b: Vector): Vector;
begin
  CheckSameLength(a, b);
  for var i := 0 to a.Length - 1 do
    a.fdata[i] += b.fdata[i];
  Result := a;
end;

static function Vector.operator -=(a, b: Vector): Vector;
begin
  CheckSameLength(a, b);
  for var i := 0 to a.Length - 1 do
    a.fdata[i] -= b.fdata[i];
  Result := a;
end;

static function Vector.operator *=(a: Vector; alpha: real): Vector;
begin
  for var i := 0 to a.Length - 1 do
    a.fdata[i] *= alpha;
  Result := a;
end;

static function Vector.operator +(v: Vector; c: real): Vector;
begin
  Result := new Vector(v.Length);
  for var i := 0 to v.Length - 1 do
    Result[i] := v[i] + c;
end;

static function Vector.operator +(c: real; v: Vector): Vector;
begin
  Result := v + c;
end;

function Vector.Sqrt: Vector := Apply(PABCSystem.Sqrt);

function Vector.Sum: real;
begin
  var s := 0.0;
  for var i := 0 to Length - 1 do
    s += fdata[i];
  Result := s;
end;

// Алиас Mean
function Vector.Average: real;
begin
  Result := Mean;
end;

function Vector.Mean: real;
begin
  CheckNonEmpty(Self);
  Result := Sum / Length;
end;

function Vector.Norm2: real;
begin
  Result := Self.Dot(Self);
end;

function Vector.Norm: real := PABCSystem.Sqrt(Norm2);

function Vector.Max: real := fdata.Max;

function Vector.Min: real := fdata.Min;

function Vector.Dot(b: Vector): real;
begin
  if Length <> b.Length then
    DimensionError(ER_DIM_MISMATCH, Length, b.Length);

  var s := 0.0;
  for var i := 0 to Length - 1 do
    s += self[i] * b[i];
  Result := s;
end;

function Vector.Apply(f: real -> real): Vector;
begin
  Result := new Vector(Length);
  for var i := 0 to Length - 1 do
    Result[i] := f(self[i]);
end;

function Vector.SubvectorBy(indices: array of integer): Vector;
begin
  var n := indices.Length;
  var resultVec := new Vector(n);

  for var i := 0 to n - 1 do
    resultVec[i] := self[indices[i]];

  Result := resultVec;
end;


//-----------------------------
//           Matrix
//-----------------------------
constructor Matrix.Create(r, c: integer);
begin
  if (r < 0) or (c < 0) then
    ArgumentOutOfRangeError(ER_MATRIX_SIZE_NEGATIVE);
  
  fdata := new real[r, c];
end;

constructor Matrix.Create(values: array[,] of real);
begin
  if values = nil then
    ArgumentNullError(ER_VALUES_NULL);
  fdata := Copy(values);
end;

constructor Matrix.Create(values: array of Vector);
begin
  if values = nil then
    ArgumentNullError(ER_VALUES_NULL);

  var r := values.Length;
  if r = 0 then
    ArgumentError(ER_VALUES_EMPTY);

  if values[0] = nil then
    ArgumentNullError(ER_VECTOR_NULL);

  var c := values[0].Length;
  if c = 0 then
    ArgumentError(ER_VECTOR_EMPTY);

  for var i := 1 to r - 1 do
  begin
    if values[i] = nil then
      ArgumentNullError(ER_VECTOR_NULL);

    if values[i].Length <> c then
      ArgumentError(ER_MATRIX_ROWS_HAVE_DIFFERENT_LENGTHS);
  end;

  fdata := new real[r, c];

  for var i := 0 to r - 1 do
    for var j := 0 to c - 1 do
      fdata[i, j] := values[i][j];
end;

function Matrix.ToArray2D: array[,] of real;
begin
  Result := Copy(fdata);
end;

function Matrix.RowToArray(r: integer): array of real;
begin
  Result := fdata.Row(r);
end;

function Matrix.Clone: Matrix;
begin
  Result := new Matrix(fdata);
end;

function Matrix.ColumnSums: Vector;
begin
  var n := RowCount;
  var p := ColCount;
  Result := new Vector(p);

  for var j := 0 to p - 1 do
    for var i := 0 to n - 1 do
      Result[j] += fdata[i, j];
end;

function Matrix.ColumnMeans: Vector;
begin
  var n := RowCount;

  if n = 0 then
    Error(ER_EMPTY_MATRIX);

  Result := ColumnSums / n;
end;

function Matrix.RowSums: Vector;
begin
  var n := RowCount;
  var p := ColCount;
  Result := new Vector(n);

  for var i := 0 to n - 1 do
    for var j := 0 to p - 1 do
      Result[i] += fdata[i, j];
end;

function Matrix.RowMeans: Vector;
begin
  var n := RowCount;
  var p := ColCount;

  if n = 0 then
    Error(ER_EMPTY_MATRIX);

  if p = 0 then
    Error(ER_EMPTY_MATRIX);

  Result := RowSums / p;
end;

function Matrix.ColumnVariances: Vector;
begin
  var n := RowCount;

  if n = 0 then
    Error(ER_EMPTY_MATRIX);

  var means := ColumnMeans;
  var p := ColCount;

  Result := new Vector(p);

  for var j := 0 to p - 1 do
    for var i := 0 to n - 1 do
    begin
      var d := fdata[i, j] - means[j];
      Result[j] += d * d;
    end;

  Result := Result / n;
end;

function Matrix.ColumnStd: Vector;
begin
  Result := ColumnVariances.Sqrt;
end;

function Matrix.RowVariances: Vector;
begin
  var n := RowCount;
  var p := ColCount;

  if n = 0 then
    Error(ER_EMPTY_MATRIX);

  if p = 0 then
    Error(ER_EMPTY_MATRIX);

  var means := RowMeans;

  Result := new Vector(n);

  for var i := 0 to n - 1 do
    for var j := 0 to p - 1 do
    begin
      var d := fdata[i, j] - means[i];
      Result[i] += d * d;
    end;

  Result := Result / p;
end;

function Matrix.RowStd: Vector;
begin
  Result := RowVariances.Sqrt;
end;

function Matrix.ColumnMins: Vector;
begin
  var n := RowCount;
  var p := ColCount;

  if n = 0 then
    Error(ER_EMPTY_MATRIX);
  
  Result := new Vector(p);

  for var j := 0 to p - 1 do
  begin
    Result[j] := fdata[0, j];

    for var i := 1 to n - 1 do
      if fdata[i, j] < Result[j] then
        Result[j] := fdata[i, j];
  end;
end;

function Matrix.ColumnMaxs: Vector;
begin
  var n := RowCount;
  var p := ColCount;

  if n = 0 then
    Error(ER_EMPTY_MATRIX);
  
  Result := new Vector(p);

  for var j := 0 to p - 1 do
  begin
    Result[j] := fdata[0, j];

    for var i := 1 to n - 1 do
      if fdata[i, j] > Result[j] then
        Result[j] := fdata[i, j];
  end;
end;

function Matrix.RowMins: Vector;
begin
  var n := RowCount;
  var p := ColCount;
  
  if n = 0 then
    Error(ER_EMPTY_MATRIX);
  
  if p = 0 then
    Error(ER_EMPTY_MATRIX);
  
  Result := new Vector(n);

  for var i := 0 to n - 1 do
  begin
    Result[i] := fdata[i, 0];

    for var j := 1 to p - 1 do
      if fdata[i, j] < Result[i] then
        Result[i] := fdata[i, j];
  end;
end;

function Matrix.RowMaxs: Vector;
begin
  var n := RowCount;
  var p := ColCount;
  
  if n = 0 then
    Error(ER_EMPTY_MATRIX);

  if p = 0 then
    Error(ER_EMPTY_MATRIX);
  
  Result := new Vector(n);

  for var i := 0 to n - 1 do
  begin
    Result[i] := fdata[i, 0];

    for var j := 1 to p - 1 do
      if fdata[i, j] > Result[i] then
        Result[i] := fdata[i, j];
  end;
end;

function Matrix.RowArgMin(i: integer): integer;
begin
  if RowCount = 0 then
    Error(ER_EMPTY_MATRIX);

  if ColCount = 0 then
    Error(ER_EMPTY_MATRIX);

  if (i < 0) or (i >= RowCount) then
    ArgumentOutOfRangeError(ER_ROW_INDEX_OUT_OF_RANGE, i, RowCount);

  var minVal := fdata[i,0];
  var arg := 0;

  for var j := 1 to ColCount - 1 do
    if fdata[i,j] < minVal then
    begin
      minVal := fdata[i,j];
      arg := j;
    end;

  Result := arg;
end;

function Matrix.RowMin(i: integer): real;
begin
  Result := fdata[i, RowArgMin(i)];
end;

function Matrix.RowArgMax(i: integer): integer;
begin
  if RowCount = 0 then
    Error(ER_EMPTY_MATRIX);

  if ColCount = 0 then
    Error(ER_EMPTY_MATRIX);

  if (i < 0) or (i >= RowCount) then
    ArgumentOutOfRangeError(ER_ROW_INDEX_OUT_OF_RANGE, i, RowCount);

  var maxVal := fdata[i,0];
  var arg := 0;

  for var j := 1 to ColCount - 1 do
    if fdata[i,j] > maxVal then
    begin
      maxVal := fdata[i,j];
      arg := j;
    end;

  Result := arg;
end;

function Matrix.RowMax(i: integer): real;
begin
  Result := fdata[i, RowArgMax(i)];
end;

function Matrix.RowSum(i: integer): real;
begin
  if RowCount = 0 then
    Error(ER_EMPTY_MATRIX);

  if (i < 0) or (i >= RowCount) then
    ArgumentOutOfRangeError(ER_ROW_INDEX_OUT_OF_RANGE, i, RowCount);

  var sum := 0.0;

  for var j := 0 to ColCount - 1 do
    sum += fdata[i,j];

  Result := sum;
end;

function Matrix.RowMean(i: integer): real;
begin
  Result := RowSum(i) / ColCount;
end;

function Matrix.RowVariance(i: integer): real;
begin
  var mean := RowMean(i);
  var sum := 0.0;

  for var j := 0 to ColCount - 1 do
  begin
    var d := fdata[i,j] - mean;
    sum += d * d;
  end;

  Result := sum / ColCount;
end;

function Matrix.RowStd(i: integer): real;
begin
  Result := Sqrt(RowVariance(i));
end;

function Matrix.ColumnArgMin(j: integer): integer;
begin
  if (j < 0) or (j >= ColCount) then
    ArgumentOutOfRangeError(ER_COL_INDEX_OUT_OF_RANGE, j, ColCount); 
    
  var n := RowCount;

  if n = 0 then
    Error(ER_EMPTY_MATRIX);

  var minVal := fdata[0,j];
  var arg := 0;

  for var i := 1 to n - 1 do
    if fdata[i,j] < minVal then
    begin
      minVal := fdata[i,j];
      arg := i;
    end;

  Result := arg;
end;

function Matrix.ColumnMin(j: integer): real;
begin
  Result := fdata[ColumnArgMin(j), j];
end;

function Matrix.ColumnArgMax(j: integer): integer;
begin
  if (j < 0) or (j >= ColCount) then
    ArgumentOutOfRangeError(ER_COL_INDEX_OUT_OF_RANGE, j, ColCount); 
    
  var n := RowCount;

  if n = 0 then
    Error(ER_EMPTY_MATRIX);

  var maxVal := fdata[0,j];
  var arg := 0;

  for var i := 1 to n - 1 do
    if fdata[i,j] > maxVal then
    begin
      maxVal := fdata[i,j];
      arg := i;
    end;

  Result := arg;
end;

function Matrix.ColumnMax(j: integer): real;
begin
  Result := fdata[ColumnArgMax(j), j];
end;

function Matrix.ColumnSum(j: integer): real;
begin
  var sum := 0.0;

  for var i := 0 to RowCount - 1 do
    sum += fdata[i,j];

  Result := sum;
end;

function Matrix.ColumnMean(j: integer): real;
begin
  Result := ColumnSum(j) / RowCount;
end;

function Matrix.ColumnVariance(j: integer): real;
begin
  var mean := ColumnMean(j);
  var sum := 0.0;

  for var i := 0 to RowCount - 1 do
  begin
    var d := fdata[i,j] - mean;
    sum += d * d;
  end;

  Result := sum / RowCount;
end;

function Matrix.ColumnStd(j: integer): real;
begin
  Result := Sqrt(ColumnVariance(j));
end;


function Matrix.FrobeniusNorm: real;
begin
  if RowCount = 0 then
    Error(ER_EMPTY_MATRIX);
  
  var s := 0.0;
  for var i := 0 to RowCount - 1 do
    for var j := 0 to ColCount - 1 do
      s += fdata[i, j] * fdata[i, j];
  Result := PABCSystem.Sqrt(s);
end;

procedure Matrix.AddScaledIdentity(lambda: real);
begin
  var n := RowCount;
  
  if RowCount <> ColCount then
    ArgumentError(ER_MATRIX_NOT_SQUARE);
  
  for var i := 0 to n - 1 do
    fdata[i, i] += lambda;
end;

static procedure Matrix.CheckSameSize(A, B: Matrix);
begin
  if (A.RowCount <> B.RowCount) or (A.ColCount <> B.ColCount) then
    DimensionError(ER_MATRIX_SIZE_MISMATCH, A.RowCount, A.ColCount, B.RowCount, B.ColCount);
end;

static procedure Matrix.CheckMulSize(A, B: Matrix);
begin
  if A.ColCount <> B.RowCount then
    DimensionError(ER_MATRIX_MUL_SIZE_MISMATCH, A.RowCount, A.ColCount, B.RowCount, B.ColCount);
end;

static procedure Matrix.CheckVecSize(A: Matrix; x: Vector);
begin
  if A.ColCount <> x.Length then
    DimensionError(ER_MATRIX_VECTOR_SIZE_MISMATCH, A.RowCount, A.ColCount, x.Length);
end;


static function Matrix.operator +(A, B: Matrix): Matrix;
begin
  CheckSameSize(A, B);
  Result := new Matrix(A.RowCount, A.ColCount);
  for var i := 0 to A.RowCount - 1 do
    for var j := 0 to A.ColCount - 1 do
      Result.fdata[i, j] := A.fdata[i, j] + B.fdata[i, j];
end;

static function Matrix.operator -(A, B: Matrix): Matrix;
begin
  CheckSameSize(A, B);
  Result := new Matrix(A.RowCount, A.ColCount);
  for var i := 0 to A.RowCount - 1 do
    for var j := 0 to A.ColCount - 1 do
      Result.fdata[i, j] := A.fdata[i, j] - B.fdata[i, j];
end;

function operator *(A: Matrix; alpha: real): Matrix; extensionmethod;
begin
  Result := new Matrix(A.RowCount, A.ColCount);
  for var i := 0 to A.RowCount - 1 do
    for var j := 0 to A.ColCount - 1 do
      Result.fdata[i, j] := alpha * A.fdata[i, j];
end;

function operator *(alpha: real; A: Matrix): Matrix; extensionmethod;
begin
  Result := A * alpha;
end;

function operator *(A: Matrix; alpha: integer): Matrix; extensionmethod;
begin
  Result := A * real(alpha);
end;

function operator *(alpha: integer; A: Matrix): Matrix; extensionmethod;
begin
  Result := real(alpha) * A;
end;


static function Matrix.operator *(A: Matrix; x: Vector): Vector;
begin
  CheckVecSize(A, x);
  Result := new Vector(A.RowCount);
  for var i := 0 to A.RowCount - 1 do
  begin
    var s := 0.0;
    for var j := 0 to A.ColCount - 1 do
      s += A.fdata[i, j] * x[j];
    Result[i] := s;
  end;
end;

{static function Matrix.operator *(A, B: Matrix): Matrix;
begin
  CheckMulSize(A, B);
  Result := new Matrix(A.RowCount, B.ColCount);
  for var i := 0 to A.RowCount - 1 do
    for var k := 0 to A.ColCount - 1 do
    begin
      var aik := A.fdata[i, k];
      if aik <> 0.0 then
        for var j := 0 to B.ColCount - 1 do
          Result.fdata[i, j] += aik * B.fdata[k, j];
    end;
end;}

static function Matrix.operator *(A, B: Matrix): Matrix;
begin
  CheckMulSize(A, B);

  var m := A.RowCount;
  var p := A.ColCount;
  var n := B.ColCount;

  var BT := B.Transpose;

  Result := new Matrix(m, n);

  for var i := 0 to m - 1 do
    for var j := 0 to n - 1 do
    begin
      var sum := 0.0;

      for var k := 0 to p - 1 do
        sum += A.fdata[i, k] * BT.fdata[j, k];

      Result.fdata[i, j] := sum;
    end;
end;

static function Matrix.operator +=(A, B: Matrix): Matrix;
begin
  CheckSameSize(A, B);
  for var i := 0 to A.RowCount - 1 do
    for var j := 0 to A.ColCount - 1 do
      A.fdata[i, j] += B.fdata[i, j];
  Result := A;
end;

static function Matrix.operator -=(A, B: Matrix): Matrix;
begin
  CheckSameSize(A, B);
  for var i := 0 to A.RowCount - 1 do
    for var j := 0 to A.ColCount - 1 do
      A.fdata[i, j] -= B.fdata[i, j];
  Result := A;
end;

static function Matrix.operator *=(A: Matrix; alpha: real): Matrix;
begin
  for var i := 0 to A.RowCount - 1 do
    for var j := 0 to A.ColCount - 1 do
      A.fdata[i, j] *= alpha;
  Result := A;
end;

function Matrix.Transpose: Matrix;
begin
  Result := new Matrix(ColCount, RowCount);
  for var i := 0 to RowCount - 1 do
    for var j := 0 to ColCount - 1 do
      Result.fdata[j, i] := fdata[i, j];
end;

function Matrix.GetRow(i: integer): Vector;
begin
  if (i < 0) or (i >= RowCount) then
    ArgumentOutOfRangeError(ER_ROW_INDEX_OUT_OF_RANGE, i, RowCount);
  
  Result := new Vector(ColCount);
  for var j := 0 to ColCount - 1 do
    Result[j] := fdata[i, j];
end;

function Matrix.TakeRows(indices: array of integer): Matrix;
begin
  if indices = nil then
    ArgumentNullError(ER_ARG_NULL, 'indices');
  
  var n := indices.Length;
  var p := ColCount;
  
  if RowCount = 0 then
    Error(ER_EMPTY_MATRIX);

  var res := new Matrix(n, p);

  var src := Data;
  var dst := res.Data;

  // --- fast path: полный срез [0..n-1]
  var isFull := n = RowCount;

  if isFull then
  begin
    for var i := 0 to n - 1 do
      if indices[i] <> i then
      begin
        isFull := False;
        break;
      end;
  end;

  if isFull then
  begin
    System.Array.Copy(src, 0, dst, 0, n * p);
    Result := res;
    exit;
  end;

  // --- обычный block-copy
  var i := 0;
  var dstOffset := 0;

  while i < n do
  begin
    var start := indices[i];
    var len := 1;

    while (i + len < n) and (indices[i + len] = start + len) do
      len += 1;

    System.Array.Copy(
      src, start * p,
      dst, dstOffset,
      len * p
    );

    dstOffset += len * p;
    i += len;
  end;

  Result := res;
end;

function Matrix.GetCol(j: integer): Vector;
begin
  if (j < 0) or (j >= ColCount) then
    ArgumentOutOfRangeError(ER_COL_INDEX_OUT_OF_RANGE, j, ColCount);
  
  Result := new Vector(RowCount);
  for var i := 0 to RowCount - 1 do
    Result[i] := fdata[i, j];
end;

static function Matrix.Identity(n: integer): Matrix;
begin
  if n < 0 then
    ArgumentOutOfRangeError(ER_MATRIX_SIZE_NEGATIVE);
  
  Result := new Matrix(n, n);
  for var i := 0 to n - 1 do
    Result[i, i] := 1.0;
end;

static function Matrix.OuterProduct(a, b: Vector): Matrix;
begin
  if a = nil then
    ArgumentNullError(ER_ARG_NULL, 'a');
  
  if b = nil then
    ArgumentNullError(ER_ARG_NULL, 'b');

  var m := a.Length;
  var n := b.Length;
  Result := new Matrix(m, n);
  
  for var i := 0 to m - 1 do
    for var j := 0 to n - 1 do
      Result[i, j] := a[i] * b[j];
end;

function Matrix.IsSymmetric(tol: real): boolean;
begin
  if RowCount <> ColCount then
  begin
    Result := false;
    exit;
  end;
  
  for var i := 0 to RowCount - 1 do
    for var j := i + 1 to ColCount - 1 do
      if Abs(fdata[i, j] - fdata[j, i]) > tol then
      begin
        Result := false;
        exit;
      end;
  
  Result := true;
end;

function Matrix.EigenSymmetric(tol: real; maxIter: integer): (Vector, Matrix);
begin
  if RowCount <> ColCount then
    ArgumentError(ER_MATRIX_NOT_SQUARE);
  
  if not IsSymmetric(tol) then
    ArgumentError(ER_MATRIX_NOT_SYMMETRIC);
  
  var n := RowCount;
  
  var M := Clone;
  var V := Matrix.Identity(n);
  
  var off := 0.0;
  for var iter := 0 to maxIter - 1 do
  begin
    // --- Норма вне-диагонали
    off := 0.0;
    for var i := 0 to n - 1 do
      for var j := i + 1 to n - 1 do
        off += M[i, j] * M[i, j];
    
    if Sqrt(off) < tol then
      break;
    
    // --- Поиск максимального вне-диагонального
    var p := 0;
    var q := 1;
    var maxVal := 0.0;
    
    for var i := 0 to n - 1 do
      for var j := i + 1 to n - 1 do
      begin
        var v1 := Abs(M[i, j]);
        if v1 > maxVal then
        begin
          maxVal := v1;
          p := i;
          q := j;
        end;
      end;
    
    var app := M[p, p];
    var aqq := M[q, q];
    var apq := M[p, q];
    
    if Abs(apq) < tol then
      continue;
    
    // --- Устойчивая формула вращения
    var tau := (aqq - app) / (2.0 * apq);
    
    var t: real;
    if tau >= 0.0 then
      t := 1.0 / (tau + Sqrt(1.0 + tau * tau))
    else
      t := -1.0 / (-tau + Sqrt(1.0 + tau * tau));
    
    var c := 1.0 / Sqrt(1.0 + t * t);
    var s := t * c;
    
    // --- Обновление M
    for var k := 0 to n - 1 do
    begin
      if (k <> p) and (k <> q) then
      begin
        var mkp := M[k, p];
        var mkq := M[k, q];
        
        var newkp := c * mkp - s * mkq;
        var newkq := s * mkp + c * mkq;
        
        M[k, p] := newkp;
        M[p, k] := newkp;
        
        M[k, q] := newkq;
        M[q, k] := newkq;
      end;
    end;
    
    M[p, p] := app - t * apq;
    M[q, q] := aqq + t * apq;
    M[p, q] := 0.0;
    M[q, p] := 0.0;
    
    // --- Обновление V
    for var k := 0 to n - 1 do
    begin
      var vkp := V[k, p];
      var vkq := V[k, q];
      
      V[k, p] := c * vkp - s * vkq;
      V[k, q] := s * vkp + c * vkq;
    end;
  end;
  
  var offNorm := Sqrt(off);

  if offNorm >= tol then
    Error(ER_EIGEN_NOT_CONVERGED, maxIter, offNorm, tol);

  
  // --- Eigenvalues
  var values := new Vector(n);
  for var i := 0 to n - 1 do
    values[i] := M[i, i];
  
  // --- Сортировка по убыванию
  for var i := 0 to n - 2 do
    for var j := i + 1 to n - 1 do
      if values[j] > values[i] then
      begin
        var tmp := values[i];
        values[i] := values[j];
        values[j] := tmp;
        
        for var k := 0 to n - 1 do
        begin
          var tv := V[k, i];
          V[k, i] := V[k, j];
          V[k, j] := tv;
        end;
      end;
  
  // --- Нормализация знаков (для стабильности PCA)
  for var i := 0 to n - 1 do
  begin
    var maxAbs := 0.0;
    var idx := 0;
    
    for var k := 0 to n - 1 do
    begin
      var v1 := Abs(V[k, i]);
      if v1 > maxAbs then
      begin
        maxAbs := v1;
        idx := k;
      end;
    end;
    
    if V[idx, i] < 0.0 then
      for var k := 0 to n - 1 do
        V[k, i] := -V[k, i];
  end;
  
  Result := (values, V);
end;

function Matrix.PCA(k: integer): (Matrix, Vector);
begin
  var m := RowCount;
  var n := ColCount;

  if k < 1 then
    ArgumentError(ER_PCA_K_INVALID);

  if k > n then
    ArgumentError(ER_PCA_K_TOO_LARGE);

  if m < 2 then
    ArgumentError(ER_PCA_NEED_TWO_SAMPLES);

  // --- Центрирование
  var Xc := Clone;

  for var j := 0 to n - 1 do
  begin
    var mean := 0.0;
    for var i := 0 to m - 1 do
      mean += Xc[i, j];
    mean /= m;

    for var i := 0 to m - 1 do
      Xc[i, j] -= mean;
  end;

  // --- Ковариационная матрица
  var C := new Matrix(n, n);

  for var i := 0 to n - 1 do
    for var j := i to n - 1 do
    begin
      var s := 0.0;
      for var t := 0 to m - 1 do
        s += Xc[t, i] * Xc[t, j];

      s /= (m - 1);

      C[i, j] := s;
      C[j, i] := s;
    end;

  // --- Eigen
  var adaptiveIter := Max(100, n * n * 10);
  var (values, V) := C.EigenSymmetric(1e-12, adaptiveIter);

  // --- Выбор k компонент
  var components := new Matrix(n, k);
  var variances := new Vector(k);

  for var j := 0 to k - 1 do
  begin
    variances[j] := values[j];
    for var i := 0 to n - 1 do
      components[i, j] := V[i, j];
  end;

  Result := (components, variances);
end;


// Helper
function Cholesky(A: Matrix): Matrix;
begin
  if A.RowCount <> A.ColCount then
    ArgumentError(ER_CHOLESKY_NOT_SQUARE);
  
  var n := A.RowCount;
  var L := new Matrix(n, n);
  
  for var i := 0 to n - 1 do
    for var j := 0 to i do
    begin
      var s := A[i, j];
      for var k := 0 to j - 1 do
        s -= L[i, k] * L[j, k];
      
      if i = j then
      begin
        if s <= 0.0 then
          raise new MLNotSPDException('Матрица не является положительно определённой (SPD)');
        L[i, i] := Sqrt(s);
      end
      else
        L[i, j] := s / L[j, j];
    end;
  
  Result := L;
end;

// Helper
function LUDecompose(A: Matrix): (Matrix, array of integer);
begin
  if A.RowCount <> A.ColCount then
    ArgumentError(ER_MATRIX_NOT_SQUARE);
  
  var n := A.RowCount;
  var LU := A.Clone;
  var p := new integer[n];
  
  for var i := 0 to n - 1 do
    p[i] := i;
  
  for var k := 0 to n - 1 do
  begin
    var maxRow := k;
    var maxVal := Abs(LU[k, k]);
    
    for var i := k + 1 to n - 1 do
      if Abs(LU[i, k]) > maxVal then
      begin
        maxVal := Abs(LU[i, k]);
        maxRow := i;
      end;
    
    if maxVal < 1e-12 then
      Error(ER_MATRIX_SINGULAR);
    
    if maxRow <> k then
    begin
      for var j := 0 to n - 1 do
      begin
        var t := LU[k, j];
        LU[k, j] := LU[maxRow, j];
        LU[maxRow, j] := t;
      end;
      
      var tp := p[k];
      p[k] := p[maxRow];
      p[maxRow] := tp;
    end;
    
    for var i := k + 1 to n - 1 do
    begin
      LU[i, k] /= LU[k, k];
      for var j := k + 1 to n - 1 do
        LU[i, j] -= LU[i, k] * LU[k, j];
    end;
  end;
  
  Result := (LU, p);
end;

// Helper
function Permute(b: Vector; p: array of integer): Vector;
begin
  var n := b.Length;
  var r := new Vector(n);
  
  for var i := 0 to n - 1 do
    r[i] := b[p[i]];
  
  Result := r;
end;

// Helper
function SolveLowerTriangularUnit(LU: Matrix; b: Vector): Vector;
begin
  var n := LU.RowCount;
  var y := new Vector(n);
  
  for var i := 0 to n - 1 do
  begin
    var s := b[i];
    for var j := 0 to i - 1 do
      s -= LU[i, j] * y[j];
    y[i] := s; 
  end;
  
  Result := y;
end;

function SolveLowerTriangular(L: Matrix; b: Vector): Vector;
begin
  var n := L.RowCount;
  var y := new Vector(n);
  
  for var i := 0 to n - 1 do
  begin
    var s := b[i];
    for var j := 0 to i - 1 do
      s -= L[i, j] * y[j];
    y[i] := s / L[i, i];
  end;
  
  Result := y;
end;


// Helper
function SolveUpperTriangular(U: Matrix; y: Vector): Vector;
begin
  var n := U.RowCount;
  var x := new Vector(n);
  
  for var i := n - 1 downto 0 do
  begin
    var s := y[i];
    for var j := i + 1 to n - 1 do
      s -= U[i, j] * x[j];
    x[i] := s / U[i, i];
  end;
  
  Result := x;
end;

function Solve(A: Matrix; b: Vector): Vector;
begin
  if A.RowCount <> A.ColCount then
    ArgumentError(ER_MATRIX_NOT_SQUARE);
  
  if b.Length <> A.RowCount then
    DimensionError(ER_VECTOR_SIZE_MISMATCH, b.Length, A.RowCount);
  
  var (LU, p) := LUDecompose(A);
  var pb := Permute(b, p);
  var y := SolveLowerTriangularUnit(LU, pb);
  Result := SolveUpperTriangular(LU, y);
end;

function SolveSPD(A: Matrix; b: Vector): Vector;
begin
  if A.RowCount <> A.ColCount then
    ArgumentError(ER_MATRIX_NOT_SQUARE);
  
  if b.Length <> A.RowCount then
    DimensionError(ER_VECTOR_SIZE_MISMATCH, b.Length, A.RowCount);
  
  var L := Cholesky(A);

  var y := SolveLowerTriangular(L, b);
  Result := SolveUpperTriangular(L.Transpose, y);
end;

function SolveAuto(A: Matrix; b: Vector): Vector;
begin
  try
    Result := SolveSPD(A, b);
  except
    on e: MLNotSPDException do
      Result := Solve(A, b);
  end;
end;

function SolveRidge(A: Matrix; b: Vector; lambda: real): Vector;
begin
  if lambda < 0.0 then
    ArgumentError(ER_LAMBDA_NEGATIVE);
  
  var m := A.RowCount;
  var n := A.ColCount;
  
  if b.Length <> m then
    DimensionError(ER_VECTOR_SIZE_MISMATCH, b.Length, m);
  
  // ------------------------------------------------------------
  // 1. Compute AtA = A^T * A
  // ------------------------------------------------------------
  var AtA := new Matrix(n, n);
  
  for var i := 0 to n - 1 do
    for var j := 0 to n - 1 do
    begin
      var s := 0.0;
      for var k := 0 to m - 1 do
        s += A[k, i] * A[k, j];
      AtA[i, j] := s;
    end;
  
  // ------------------------------------------------------------
  // 2. Add lambda * I
  // ------------------------------------------------------------
  if lambda <> 0.0 then
    for var i := 0 to n - 1 do
      AtA[i, i] += lambda;
  
  // ------------------------------------------------------------
  // 3. Compute Atb = A^T * b
  // ------------------------------------------------------------
  var Atb := new Vector(n);
  
  for var i := 0 to n - 1 do
  begin
    var s := 0.0;
    for var k := 0 to m - 1 do
      s += A[k, i] * b[k];
    Atb[i] := s;
  end;
  
  // ------------------------------------------------------------
  // 4. Solve SPD system
  // ------------------------------------------------------------
  Result := SolveSPD(AtA, Atb);
end;


function SolveLeastSquaresQR(A: Matrix; b: Vector): Vector;
begin
  if A.RowCount <> b.Length then
    DimensionError(ER_DIM_MISMATCH, A.RowCount, b.Length);

  var m := A.RowCount;
  var n := A.ColCount;

  if m < n then
    ArgumentError(ER_QR_REQUIRES_M_GE_N);

  var R := A.Clone;
  var y := b.Clone;

  for var k := 0 to n - 1 do
  begin
    // ---- вычислить норму столбца k начиная с строки k
    var normx := 0.0;
    for var i := k to m - 1 do
      normx += R[i,k] * R[i,k];
    
    normx := Sqrt(normx);
    
    if normx = 0 then
      continue;
    
    if R[k,k] >= 0 then
      normx := -normx;
    
    R[k,k] -= normx;
    
    var vnorm2 := 0.0;
    for var i := k to m - 1 do
      vnorm2 += R[i,k] * R[i,k];
    
    if vnorm2 < 1e-12 then
      continue;
    
    var beta := 2.0 / vnorm2;

    // ---- применить отражение к R
    for var j := k + 1 to n - 1 do
    begin
      var s := 0.0;
      for var i := k to m - 1 do
        s += R[i,k] * R[i,j];

      s *= beta;

      for var i := k to m - 1 do
        R[i,j] -= s * R[i,k];
    end;

    // ---- применить отражение к y
    var sy := 0.0;
    for var i := k to m - 1 do
      sy += R[i,k] * y[i];

    sy *= beta;

    for var i := k to m - 1 do
      y[i] -= sy * R[i,k];

    // ---- восстановить диагональный элемент
    R[k,k] := normx;

    for var i := k + 1 to m - 1 do
      R[i,k] := 0.0;
  end;

  // ---- Back substitution (решаем R[0:n,0:n] * x = y[0:n])
  var x := new Vector(n);

  for var i := n - 1 downto 0 do
  begin
    var sum := y[i];

    for var j := i + 1 to n - 1 do
      sum -= R[i,j] * x[j];

    if Abs(R[i,i]) < 1e-14 then
      ArgumentError(ER_SINGULAR_MATRIX);

    x[i] := sum / R[i,i];
  end;

  Result := x;
end;


end.