type ConstructorArg =
  | number[][]
  | number[]
  | number
  | Float64Array
  | Matrix
  | MatrixOptions;
type MatrixLike = Matrix | number;

interface MatrixOptions {
  array?: Float64Array;
  rows: number;
  cols: number;
  offset?: number;
  rowStride?: number;
  colStride?: number;
}

function is2DArray<T>(array: T[] | T[][]): array is T[][] {
  return Array.isArray(array[0]);
}

/**
 * A custom-made high-performance two-dimensional matrix backed by a
 * `Float64Array`.
 *
 * Most operations are not checked for correctness. The consumer of this class
 * is expected to ensure correctness of all operations. Invalid operations
 * (e.g. out-of-bounds array access, matrix multiplication with incompatible
 * sizes) will result in undefined behavior.
 *
 * @author Rishov Sarkar
 */
export class Matrix {
  private _array: Float64Array;
  private _rows: number;
  private _cols: number;
  private _offset: number;
  private _rowStride: number;
  private _colStride: number;

  /**
   * The number of rows in this matrix.
   * @see `this.reshape` to change the number of rows and columns.
   */
  get rows() {
    return this._rows;
  }

  /**
   * The number of columns in this matrix.
   * @see `this.reshape` to change the number of rows and columns.
   */
  get cols() {
    return this._cols;
  }

  /**
   * Create a 0&times;0 matrix.
   */
  constructor();
  /**
   * Create a `size`&times;`size` matrix of all zeros.
   */
  constructor(size: number);
  /**
   * Create a 1&times;`row.length` row vector populated from the given `row`.
   */
  constructor(row: number[]);
  /**
   * Create an `array.length`&times;`array[0].length` matrix populated from
   * the given 2-D `array`.
   */
  constructor(array: number[][]);
  /**
   * Create a new `Matrix` from the given `Matrix`.
   *
   * Note that the new `Matrix` shares the same backing array as the provided
   * `Matrix`.
   */
  constructor(matrix: Matrix);
  /**
   * Create a 1&times;`row.length` row vector using the given `row` as the
   * backing array.
   *
   * Note that the given `row` is used directly as a backing array; its
   * contents are not copied.
   */
  constructor(row: Float64Array);
  /**
   * Create a new `Matrix` with the options specified.
   *
   * This constructor form is the only one that allows full customization of
   * all the `Matrix` internals, including offset and strides.
   */
  constructor(options: MatrixOptions);
  /**
   * Create a `rows`&times;`cols` matrix of all zeros.
   */
  constructor(rows: number, cols: number);
  /**
   * Create a `rows`&times;`cols` matrix using the given `array` as the backing
   * array.
   *
   * Note that the given `array` is used directly as a backing array; its
   * contents are not copied.
   */
  constructor(rows: number, cols: number, array: Float64Array);
  constructor(...args: ConstructorArg[]) {
    switch (args.length) {
      case 0: {
        this._rows = 0;
        this._cols = 0;
        this._array = new Float64Array();
        this._offset = 0;
        this._rowStride = 0;
        this._colStride = 0;
        break;
      }
      case 1: {
        const [arg] = args;
        if (typeof arg === 'number') {
          this._rows = arg;
          this._cols = arg;
          this._array = new Float64Array(arg * arg);
          this._offset = 0;
          this._rowStride = this._cols;
          this._colStride = 1;
        } else if (Array.isArray(arg)) {
          if (is2DArray<number>(arg)) {
            this._rows = arg.length;
            this._cols = arg[0].length;
            this._array = new Float64Array(this._rows * this._cols);
            let index = 0;
            for (let i = 0; i < this._rows; i++) {
              for (let j = 0; j < this._cols; j++) {
                this._array[index] = arg[i][j];
                index++;
              }
            }
          } else {
            this._rows = 1;
            this._cols = arg.length;
            this._array = new Float64Array(arg);
          }
          this._offset = 0;
          this._rowStride = this._cols;
          this._colStride = 1;
        } else if (arg instanceof Matrix) {
          this._rows = arg._rows;
          this._cols = arg._cols;
          this._array = arg._array;
          this._offset = arg._offset;
          this._rowStride = arg._rowStride;
          this._colStride = arg._colStride;
        } else if (arg instanceof Float64Array) {
          this._rows = 1;
          this._cols = arg.length;
          this._array = arg;
          this._offset = 0;
          this._rowStride = this._cols;
          this._colStride = 1;
        } else {
          this._rows = arg.rows;
          this._cols = arg.cols;
          this._array = arg.array ?? new Float64Array(this._rows * this._cols);
          this._offset = arg.offset ?? 0;
          this._rowStride = arg.rowStride ?? 0;
          this._colStride = arg.colStride ?? 1;
        }
        break;
      }
      case 2: {
        const [rows, cols] = args;
        if (typeof rows !== 'number') {
          throw new Error('rows must be a number');
        }
        if (typeof cols !== 'number') {
          throw new Error('cols must be a number');
        }
        this._rows = rows;
        this._cols = cols;
        this._array = new Float64Array(rows * cols);
        this._offset = 0;
        this._rowStride = this._cols;
        this._colStride = 1;
        break;
      }
      case 3: {
        const [rows, cols, array] = args;
        if (typeof rows !== 'number') {
          throw new Error('rows must be a number');
        }
        if (typeof cols !== 'number') {
          throw new Error('cols must be a number');
        }
        if (!(array instanceof Float64Array)) {
          throw new Error('array must be a Float64Array');
        }
        this._rows = rows;
        this._cols = cols;
        this._array = array;
        this._offset = 0;
        this._rowStride = this._cols;
        this._colStride = 1;
        break;
      }
      default:
        throw new Error('no constructor matching number of arguments');
    }
  }

  /**
   * Deeply copy this matrix. This matrix is not modified.
   * @returns a deep copy of this matrix.
   */
  copy() {
    const copy = new Matrix(this);
    copy._array = new Float64Array(copy._array);
    return copy;
  }

  private _index(row: number, col?: number) {
    if (col === undefined) {
      const stride = this._rows !== 1 ? this._rowStride : this._colStride;
      return row * stride + this._offset;
    }
    return row * this._rowStride + col * this._colStride + this._offset;
  }

  /**
   * Get the value at the specified `index` of a row or column vector.
   * @param index column to access if this is a row vector, or row to access if
   * this is a column vector
   */
  at(index: number): number;
  /**
   * Get the value at the specified `row` and `col`.
   */
  at(row: number, col: number): number;
  at(row: number, col?: number) {
    return this._array[this._index(row, col)];
  }

  /**
   * Set the value at the specified `index` of a row or column vector to
   * `value`.
   * @param index column to access if this is a row vector, or row to access if
   * this is a column vector
   */
  set(index: number, value: number): void;
  /**
   * Set the value at the specified `row` and `col` to `value`.
   */
  set(row: number, col: number, value: number): void;
  set(row: number, colOrValue: number, value?: number) {
    const col = value !== undefined ? colOrValue : undefined;
    this._array[this._index(row, col)] = value ?? colOrValue;
  }

  /**
   * Get one row of this matrix. The returned row is backed by the same backing
   * array as this matrix. This is a constant-time, constant-space operation.
   * @returns a new `Matrix` containing only the specified `row`
   */
  row(row: number) {
    const result = new Matrix(this);
    result._rows = 1;
    result._offset += row * result._rowStride;
    return result;
  }

  /**
   * Get an array of all rows of this matrix. The returned rows are backed by
   * the same backing array as this matrix. This operation takes time and space
   * linear with respect to the number of rows.
   * @returns an array of `Matrix`s containing each row
   */
  allRows() {
    const result = new Array<Matrix>(this._rows);
    for (let i = 0; i < this._rows; i++) {
      result[i] = this.row(i);
    }
    return result;
  }

  /**
   * Get one column of this matrix. The returned column is backed by the same
   * backing array as this matrix. This is a constant-time, constant-space
   * operation.
   * @returns a new `Matrix` containing only the specified `col`
   */
  col(col: number) {
    const result = new Matrix(this);
    result._cols = 1;
    result._offset += col * result._colStride;
    return result;
  }

  /**
   * Get an array of all columns of this matrix. The returned columns are
   * backed by the same backing array as this matrix. This operation takes time
   * and space linear with respect to the number of columns.
   * @returns an array of `Matrix`s containing each column
   */
  allCols() {
    const result = new Array<Matrix>(this._cols);
    for (let i = 0; i < this._cols; i++) {
      result[i] = this.col(i);
    }
    return result;
  }

  slice(fromRow = 0, toRow = this._rows, fromCol = 0, toCol = this._cols) {
    const result = new Matrix(this);
    result._rows = toRow - fromRow;
    result._cols = toCol - fromCol;
    result._offset = fromRow * result._rowStride + fromCol * result._colStride;
    return result;
  }

  sliceRows(fromRow?: number, toRow?: number) {
    return this.slice(fromRow, toRow, undefined, undefined);
  }

  sliceCols(fromCol?: number, toCol?: number) {
    return this.slice(undefined, undefined, fromCol, toCol);
  }

  transpose() {
    const prevRowStride = this._rowStride;
    this._rowStride = this._colStride;
    this._colStride = prevRowStride;
    const prevRows = this._rows;
    this._rows = this._cols;
    this._cols = prevRows;
    return this;
  }

  static transpose(matrix: Matrix) {
    return new Matrix(matrix).transpose();
  }

  optimize(colMajor = false) {
    const optDenseStride = 1;
    const optSparseStride = colMajor ? this._rows : this._cols;
    const denseStride = colMajor ? this._rowStride : this._colStride;
    const sparseStride = colMajor ? this._colStride : this._rowStride;
    const numDense = colMajor ? this._rows : this._cols;
    const numSparse = colMajor ? this._cols : this._rows;
    const numElements = numDense * numSparse;

    if (denseStride === optDenseStride && sparseStride === optSparseStride) {
      if (this._offset !== 0 || this._array.length !== numElements) {
        this._array = this._array.subarray(
          this._offset,
          this._offset + numElements
        );
        this._offset = 0;
      }
      return this;
    }

    const array = new Float64Array(numElements);
    let srcSparseOffset = this._offset;
    let dstOffset = 0;
    for (let i = 0; i < numSparse; i++) {
      let srcOffset = srcSparseOffset;
      for (let j = 0; j < numDense; j++) {
        array[dstOffset] = this._array[srcOffset];
        srcOffset += denseStride;
        dstOffset++;
      }
      srcSparseOffset += sparseStride;
    }

    this._array = array;
    this._offset = 0;
    this._rowStride = colMajor ? 1 : this._cols;
    this._colStride = colMajor ? this._rows : 1;
    return this;
  }

  reshape(rows?: number, cols?: number, colMajor = false) {
    const numElements = this._rows * this._cols;
    rows = rows ?? numElements / (cols ?? numElements);
    cols = cols ?? numElements / rows;
    this.optimize(colMajor);
    this._rows = rows;
    this._cols = cols;
    this._rowStride = colMajor ? 1 : cols;
    this._colStride = colMajor ? rows : 1;
    return this;
  }

  static reshape(
    matrix: Matrix,
    rows?: number,
    cols?: number,
    colMajor = false
  ) {
    return new Matrix(matrix).reshape(rows, cols, colMajor);
  }

  multiply(matrix: Matrix, optimize = false) {
    const first = optimize ? new Matrix(this).optimize(false) : this;
    const second = optimize ? new Matrix(matrix).optimize(true) : matrix;
    const array = new Float64Array(first._rows * second._cols);
    let firstRowOffset = first._offset;
    let resultOffset = 0;
    for (let i = 0; i < first._rows; i++) {
      let secondColOffset = second._offset;
      for (let j = 0; j < second._cols; j++) {
        let firstOffset = firstRowOffset;
        let secondOffset = secondColOffset;
        let dot = 0;
        for (let k = 0; k < first._cols; k++) {
          dot += first._array[firstOffset] * second._array[secondOffset];
          firstOffset += first._colStride;
          secondOffset += second._rowStride;
        }
        array[resultOffset] = dot;
        secondColOffset += second._colStride;
        resultOffset++;
      }
      firstRowOffset += first._rowStride;
    }
    this._array = array;
    this._cols = second._cols;
    this._rowStride = this._cols;
    this._colStride = 1;
    return this;
  }

  static multiply(first: Matrix, second: Matrix, optimize = false) {
    return new Matrix(first).multiply(second, optimize);
  }

  magnitude() {
    return Math.sqrt(Matrix.dot(this, this));
  }

  normalize() {
    return Matrix._normalize(this, this);
  }

  static normalize(vector: Matrix) {
    return Matrix._normalize(undefined, vector);
  }

  private static _normalize(dst: Matrix | undefined, src: MatrixLike) {
    src = Matrix._broadcast(src, 0);
    dst = !dst || !Matrix._fits(dst, src) ? Matrix.zerosLike(src) : dst;
    const mag = src.magnitude();
    const length = src._rows !== 1 ? src._rows : src._cols;
    const srcStride = src._rows !== 1 ? src._rowStride : src._colStride;
    const dstStride = dst._rows !== 1 ? dst._rowStride : dst._colStride;
    let srcOffset = src._offset;
    let dstOffset = dst._offset;
    for (let i = 0; i < length; i++) {
      dst._array[dstOffset] = src._array[srcOffset] / mag;
      srcOffset += srcStride;
      dstOffset += dstStride;
    }
    return dst;
  }

  negate() {
    return this.scale(-1);
  }

  static negate(matrix: Matrix) {
    return Matrix.scale(matrix, -1);
  }

  static dot(a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    const length = a._rows !== 1 ? a._rows : a._cols;
    const aStride = a._rows !== 1 ? a._rowStride : a._colStride;
    const bStride = b._rows !== 1 ? b._rowStride : b._colStride;
    let aOffset = a._offset;
    let bOffset = b._offset;
    let result = 0;
    for (let i = 0; i < length; i++) {
      result += a._array[aOffset] * b._array[bOffset];
      aOffset += aStride;
      bOffset += bStride;
    }
    return result;
  }

  cross(vector: MatrixLike) {
    return Matrix._cross(this, this, vector);
  }

  static cross(first: MatrixLike, second: MatrixLike) {
    return Matrix._cross(undefined, first, second);
  }

  private static _cross(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;

    const aStride = a._rows !== 1 ? a._rowStride : a._colStride;
    let aOffset = a._offset;
    const a1 = a._array[aOffset];
    aOffset += aStride;
    const a2 = a._array[aOffset];
    aOffset += aStride;
    const a3 = a._array[aOffset];

    const bStride = b._rows !== 1 ? b._rowStride : b._colStride;
    let bOffset = b._offset;
    const b1 = b._array[bOffset];
    bOffset += bStride;
    const b2 = b._array[bOffset];
    bOffset += bStride;
    const b3 = b._array[bOffset];

    const dstStride = b._rows !== 1 ? b._rowStride : b._colStride;
    let dstOffset = dst._offset;
    dst._array[dstOffset] = a2 * b3 - a3 * b2;
    dstOffset += dstStride;
    dst._array[dstOffset] = a3 * b1 - a1 * b3;
    dstOffset += dstStride;
    dst._array[dstOffset] = a1 * b2 - a2 * b1;
    return dst;
  }

  add(matrix: MatrixLike) {
    return Matrix._add(this, this, matrix);
  }

  static add(first: MatrixLike, second: MatrixLike) {
    return Matrix._add(undefined, first, second);
  }

  private static _add(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = a._array[aOffset] + b._array[bOffset];
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  subtract(matrix: MatrixLike) {
    return Matrix._subtract(this, this, matrix);
  }

  static subtract(first: MatrixLike, second: MatrixLike) {
    return Matrix._subtract(undefined, first, second);
  }

  private static _subtract(
    dst: Matrix | undefined,
    a: MatrixLike,
    b: MatrixLike
  ) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = a._array[aOffset] - b._array[bOffset];
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  scale(matrix: MatrixLike) {
    return Matrix._scale(this, this, matrix);
  }

  static scale(first: MatrixLike, second: MatrixLike) {
    return Matrix._scale(undefined, first, second);
  }

  private static _scale(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = a._array[aOffset] * b._array[bOffset];
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  divide(matrix: MatrixLike) {
    return Matrix._divide(this, this, matrix);
  }

  static divide(first: MatrixLike, second: MatrixLike) {
    return Matrix._divide(undefined, first, second);
  }

  private static _divide(
    dst: Matrix | undefined,
    a: MatrixLike,
    b: MatrixLike
  ) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = a._array[aOffset] / b._array[bOffset];
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  gt(matrix: MatrixLike) {
    return Matrix._gt(this, this, matrix);
  }

  static gt(first: MatrixLike, second: MatrixLike) {
    return Matrix._gt(undefined, first, second);
  }

  private static _gt(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = +(a._array[aOffset] > b._array[bOffset]);
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  lt(matrix: MatrixLike) {
    return Matrix._lt(this, this, matrix);
  }

  static lt(first: MatrixLike, second: MatrixLike) {
    return Matrix._lt(undefined, first, second);
  }

  private static _lt(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = +(a._array[aOffset] < b._array[bOffset]);
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  gte(matrix: MatrixLike) {
    return Matrix._gte(this, this, matrix);
  }

  static gte(first: MatrixLike, second: MatrixLike) {
    return Matrix._gte(undefined, first, second);
  }

  private static _gte(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = +(a._array[aOffset] >= b._array[bOffset]);
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  lte(matrix: MatrixLike) {
    return Matrix._lte(this, this, matrix);
  }

  static lte(first: MatrixLike, second: MatrixLike) {
    return Matrix._lte(undefined, first, second);
  }

  private static _lte(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = +(a._array[aOffset] <= b._array[bOffset]);
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  eq(matrix: MatrixLike) {
    return Matrix._eq(this, this, matrix);
  }

  static eq(first: MatrixLike, second: MatrixLike) {
    return Matrix._eq(undefined, first, second);
  }

  private static _eq(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = +(a._array[aOffset] === b._array[bOffset]);
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  mask(vector: Matrix) {
    return Matrix._mask(this, this, vector);
  }

  static mask(matrix: Matrix, vector: Matrix) {
    return Matrix._mask(undefined, matrix, vector);
  }

  private static _mask(dst: Matrix | undefined, a: Matrix, b: Matrix) {
    dst = dst ?? Matrix.zerosLike(a);
    const maskRows = b._rows !== 1;
    const numOuter = maskRows ? a._rows : a._cols;
    const numInner = maskRows ? a._cols : a._rows;
    const aOuterStride = maskRows ? a._rowStride : a._colStride;
    const aInnerStride = maskRows ? a._colStride : a._rowStride;
    const bStride = maskRows ? b._rowStride : b._colStride;
    const dstOuterStride = maskRows ? dst._rowStride : dst._colStride;
    const dstInnerStride = maskRows ? dst._colStride : dst._rowStride;

    let aOuterOffset = a._offset;
    let bOffset = b._offset;
    let dstOuterOffset = dst._offset;
    let numNotMasked = 0;
    for (let i = 0; i < numOuter; i++) {
      let aOffset = aOuterOffset;
      let dstOffset = dstOuterOffset;
      if (b._array[bOffset]) {
        for (let j = 0; j < numInner; j++) {
          dst._array[dstOffset] = a._array[aOffset];
          aOffset += aInnerStride;
          dstOffset += dstInnerStride;
        }
        dstOuterOffset += dstOuterStride;
        numNotMasked++;
      }
      aOuterOffset += aOuterStride;
      bOffset += bStride;
    }

    if (maskRows) {
      dst._rows = numNotMasked;
    } else {
      dst._cols = numNotMasked;
    }
    return dst;
  }

  pow(matrix: MatrixLike) {
    return Matrix._pow(this, this, matrix);
  }

  static pow(first: MatrixLike, second: MatrixLike) {
    return Matrix._pow(undefined, first, second);
  }

  private static _pow(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        dst._array[dstOffset] = a._array[aOffset] ** b._array[bOffset];
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  min(matrix: MatrixLike) {
    return Matrix._min(this, this, matrix);
  }

  static min(first: MatrixLike, second: MatrixLike) {
    return Matrix._min(undefined, first, second);
  }

  private static _min(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        const aValue = a._array[aOffset];
        const bValue = b._array[bOffset];
        dst._array[dstOffset] = aValue < bValue ? aValue : bValue;
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  max(matrix: MatrixLike) {
    return Matrix._max(this, this, matrix);
  }

  static max(first: MatrixLike, second: MatrixLike) {
    return Matrix._max(undefined, first, second);
  }

  private static _max(dst: Matrix | undefined, a: MatrixLike, b: MatrixLike) {
    a = Matrix._broadcast(a, b);
    b = Matrix._broadcast(b, a);
    dst = !dst || !Matrix._fits(dst, a) ? Matrix.zerosLike(a) : dst;
    let aRowOffset = a._offset;
    let bRowOffset = b._offset;
    let dstRowOffset = dst._offset;
    for (let i = 0; i < a._rows; i++) {
      let aOffset = aRowOffset;
      let bOffset = bRowOffset;
      let dstOffset = dstRowOffset;
      for (let j = 0; j < a._cols; j++) {
        const aValue = a._array[aOffset];
        const bValue = b._array[bOffset];
        dst._array[dstOffset] = aValue > bValue ? aValue : bValue;
        aOffset += a._colStride;
        bOffset += b._colStride;
        dstOffset += dst._colStride;
      }
      aRowOffset += a._rowStride;
      bRowOffset += b._rowStride;
      dstRowOffset += dst._rowStride;
    }
    return dst;
  }

  sumRows() {
    return Matrix._sumRows(this, this);
  }

  static sumRows(matrix: Matrix) {
    return Matrix._sumRows(undefined, matrix);
  }

  private static _sumRows(dst: Matrix | undefined, src: Matrix) {
    dst = dst ?? Matrix.zerosLike(src);
    let srcColOffset = src._offset;
    let dstOffset = dst._offset;
    for (let i = 0; i < src._cols; i++) {
      let srcOffset = srcColOffset;
      let sum = 0;
      for (let j = 0; j < src._rows; j++) {
        sum += src._array[srcOffset];
        srcOffset += src._rowStride;
      }
      dst._array[dstOffset] = sum;
      srcColOffset += src._colStride;
      dstOffset += dst._colStride;
    }
    dst._rows = 1;
    return dst;
  }

  sumCols() {
    return Matrix._sumCols(this, this);
  }

  static sumCols(matrix: Matrix) {
    return Matrix._sumCols(undefined, matrix);
  }

  private static _sumCols(dst: Matrix | undefined, src: Matrix) {
    dst = dst ?? Matrix.zerosLike(src);
    let srcRowOffset = src._offset;
    let dstOffset = dst._offset;
    for (let i = 0; i < src._rows; i++) {
      let srcOffset = srcRowOffset;
      let sum = 0;
      for (let j = 0; j < src._cols; j++) {
        sum += src._array[srcOffset];
        srcOffset += src._colStride;
      }
      dst._array[dstOffset] = sum;
      srcRowOffset += src._rowStride;
      dstOffset += dst._rowStride;
    }
    dst._cols = 1;
    return dst;
  }

  flipHorizontal() {
    this._offset += (this._cols - 1) * this._colStride;
    this._colStride *= -1;
    return this;
  }

  static flipHorizontal(matrix: Matrix) {
    return new Matrix(matrix).flipHorizontal();
  }

  flipVertical() {
    this._offset += (this._rows - 1) * this._rowStride;
    this._rowStride *= -1;
    return this;
  }

  static flipVertical(matrix: Matrix) {
    return new Matrix(matrix).flipVertical();
  }

  reverse() {
    return this.flipHorizontal().flipVertical();
  }

  static reverse(matrix: Matrix) {
    return new Matrix(matrix).reverse();
  }

  stackHorizontal(matrix: MatrixLike) {
    return Matrix._stack(this, this, matrix, true);
  }

  static stackHorizontal(first: MatrixLike, second: MatrixLike) {
    return Matrix._stack(undefined, first, second, true);
  }

  stackVertical(matrix: MatrixLike) {
    return Matrix._stack(this, this, matrix, false);
  }

  static stackVertical(first: MatrixLike, second: MatrixLike) {
    return Matrix._stack(undefined, first, second, false);
  }

  private static _stack(
    dst: Matrix | undefined,
    a: MatrixLike,
    b: MatrixLike,
    horizontal: boolean
  ) {
    a = Matrix._broadcast(a, b, horizontal, !horizontal);
    b = Matrix._broadcast(b, a, horizontal, !horizontal);
    const rows = horizontal ? a._rows : a._rows + b._rows;
    const cols = horizontal ? a._cols + b._cols : a._cols;
    const array = new Float64Array(rows * cols);
    dst = dst ?? new Matrix(rows, cols, array);

    const numDense = horizontal ? rows : cols;
    const aNumSparse = horizontal ? a._cols : a._rows;
    const bNumSparse = horizontal ? b._cols : b._rows;
    const aDenseStride = horizontal ? a._rowStride : a._colStride;
    const aSparseStride = horizontal ? a._colStride : a._rowStride;
    const bDenseStride = horizontal ? b._rowStride : b._colStride;
    const bSparseStride = horizontal ? b._colStride : b._rowStride;
    let dstOffset = 0;

    let aSparseOffset = a._offset;
    for (let i = 0; i < aNumSparse; i++) {
      let aOffset = aSparseOffset;
      for (let j = 0; j < numDense; j++) {
        array[dstOffset] = a._array[aOffset];
        aOffset += aDenseStride;
        dstOffset++;
      }
      aSparseOffset += aSparseStride;
    }

    let bSparseOffset = b._offset;
    for (let i = 0; i < bNumSparse; i++) {
      let bOffset = bSparseOffset;
      for (let j = 0; j < numDense; j++) {
        array[dstOffset] = b._array[bOffset];
        bOffset += bDenseStride;
        dstOffset++;
      }
      bSparseOffset += bSparseStride;
    }

    dst._array = array;
    dst._rows = rows;
    dst._cols = cols;
    dst._offset = 0;
    dst._rowStride = horizontal ? 1 : cols;
    dst._colStride = horizontal ? rows : 1;
    return dst;
  }

  toString() {
    let result = `Matrix (${this._rows}\xD7${this._cols}) [`;
    if (this._rows !== 0 && this._cols !== 0) {
      let rowOffset = this._offset;
      for (let i = 0; i < this._rows; i++) {
        result += '\n  ';
        let offset = rowOffset;
        for (let j = 0; j < this._cols; j++) {
          if (j !== 0) {
            result += ', ';
          }
          result += this._array[offset];
          offset += this._colStride;
        }
        rowOffset += this._rowStride;
      }
      result += '\n';
    }
    result += ']';
    return result;
  }

  toArray(twoDimensional: false): number[];
  toArray(twoDimensional?: true): number[][];
  toArray(twoDimensional = true) {
    return twoDimensional ? this._to2DArray() : this._to1DArray();
  }

  private _to2DArray() {
    const result = new Array<number[]>(this._rows);
    let rowOffset = this._offset;
    for (let i = 0; i < this._rows; i++) {
      const resultRow = new Array<number>(this._cols);
      let offset = rowOffset;
      for (let j = 0; j < this._cols; j++) {
        resultRow[j] = this._array[offset];
        offset += this._colStride;
      }
      result[i] = resultRow;
      rowOffset += this._rowStride;
    }
    return result;
  }

  private _to1DArray() {
    const result = new Array<number>(this._rows * this._cols);
    let rowOffset = this._offset;
    let resultIndex = 0;
    for (let i = 0; i < this._rows; i++) {
      let offset = rowOffset;
      for (let j = 0; j < this._cols; j++) {
        result[resultIndex] = this._array[offset];
        offset += this._colStride;
        resultIndex++;
      }
      rowOffset += this._rowStride;
    }
    return result;
  }

  static identity(size: number) {
    const result = new Matrix(size);
    const stride = result._rowStride + result._colStride;
    let offset = result._offset;
    for (let i = 0; i < size; i++) {
      result._array[offset] = 1;
      offset += stride;
    }
    return result;
  }

  static zerosLike(matrix: Matrix) {
    return new Matrix(matrix._rows, matrix._cols);
  }

  private static _broadcast(
    matrix: MatrixLike,
    other: MatrixLike,
    matchRows = true,
    matchCols = true
  ) {
    if (typeof matrix === 'number') {
      return new Matrix({
        rows: matchRows && typeof other !== 'number' ? other._rows : 1,
        cols: matchCols && typeof other !== 'number' ? other._cols : 1,
        offset: 0,
        rowStride: 0,
        colStride: 0,
        array: Float64Array.of(matrix),
      });
    }
    if (typeof other === 'number') {
      return matrix;
    }
    let result: Matrix | undefined;
    if (matchRows && matrix._rows === 1 && other._rows !== 1) {
      result = result ?? new Matrix(matrix);
      result._rows = other._rows;
      result._rowStride = 0;
    }
    if (matchCols && matrix._cols === 1 && other._cols !== 1) {
      result = result ?? new Matrix(matrix);
      result._cols = other._cols;
      result._colStride = 0;
    }
    return result ?? matrix;
  }

  private static _fits(dst: Matrix, src: Matrix) {
    return (
      dst._rows >= src._rows &&
      dst._cols >= src._cols &&
      (src._rows <= 1 || dst._rowStride !== 0) &&
      (src._cols <= 1 || dst._colStride !== 0)
    );
  }
}
