Class Matrix < Object

# File ../ruby/lib/matrix.rb, line 139
def Matrix.[](*rows)
  Matrix.rows(rows, false)
end

# File ../ruby/lib/matrix.rb, line 184
def Matrix.build(row_size, column_size = row_size)
  row_size = CoercionHelper.coerce_to_int(row_size)
  column_size = CoercionHelper.coerce_to_int(column_size)
  raise ArgumentError if row_size < 0 || column_size < 0
  return to_enum :build, row_size, column_size unless block_given?
  rows = Array.new(row_size) do |i|
    Array.new(column_size) do |j|
      yield i, j
    end
  end
  new rows, column_size
end

# File ../ruby/lib/matrix.rb, line 269
def Matrix.column_vector(column)
  column = convert_to_array(column)
  new [column].transpose, 1
end

# File ../ruby/lib/matrix.rb, line 169
def Matrix.columns(columns)
  Matrix.rows(columns, false).transpose
end

# File ../ruby/lib/matrix.rb, line 204
def Matrix.diagonal(*values)
  size = values.size
  rows = Array.new(size) {|j|
    row = Array.new(size, 0)
    row[j] = values[j]
    row
  }
  new rows
end

# File ../ruby/lib/matrix.rb, line 287
def Matrix.empty(row_size = 0, column_size = 0)
  Matrix.Raise ArgumentError, "One size must be 0" if column_size != 0 && row_size != 0
  Matrix.Raise ArgumentError, "Negative size" if column_size < 0 || row_size < 0

  new([[]]*row_size, column_size)
end

# File ../ruby/lib/matrix.rb, line 231
def Matrix.identity(n)
  Matrix.scalar(n, 1)
end

# File ../ruby/lib/matrix.rb, line 297
def initialize(rows, column_size = rows[0].size)
  # No checking is done at this point. rows must be an Array of Arrays.
  # column_size must be the size of the first row, if there is one,
  # otherwise it *must* be specified and can be any integer >= 0
  @rows = rows
  @column_size = column_size
end

# File ../ruby/lib/matrix.rb, line 256
def Matrix.row_vector(row)
  row = convert_to_array(row)
  new [row]
end

# File ../ruby/lib/matrix.rb, line 151
def Matrix.rows(rows, copy = true)
  rows = convert_to_array(rows)
  rows.map! do |row|
    convert_to_array(row, copy)
  end
  size = (rows[0] || []).size
  rows.each do |row|
    Matrix.Raise ErrDimensionMismatch, "row size differs (#{row.size} should be #{size})" unless row.size == size
  end
  new rows, size
end

# File ../ruby/lib/matrix.rb, line 221
def Matrix.scalar(n, value)
  Matrix.diagonal(*Array.new(n, value))
end

# File ../ruby/lib/matrix.rb, line 245
def Matrix.zero(row_size, column_size = row_size)
  rows = Array.new(row_size){Array.new(column_size, 0)}
  new rows, column_size
end

# File ../ruby/lib/matrix.rb, line 802
def *(m) # m is matrix or vector or number
  case(m)
  when Numeric
    rows = @rows.collect {|row|
      row.collect {|e| e * m }
    }
    return new_matrix rows, column_size
  when Vector
    m = Matrix.column_vector(m)
    r = self * m
    return r.column(0)
  when Matrix
    Matrix.Raise ErrDimensionMismatch if column_size != m.row_size

    rows = Array.new(row_size) {|i|
      Array.new(m.column_size) {|j|
        (0 ... column_size).inject(0) do |vij, k|
          vij + self[i, k] * m[k, j]
        end
      }
    }
    return new_matrix rows, m.column_size
  else
    return apply_through_coercion(m, __method__)
  end
end

# File ../ruby/lib/matrix.rb, line 969
def ** (other)
  case other
  when Integer
    x = self
    if other <= 0
      x = self.inverse
      return Matrix.identity(self.column_size) if other == 0
      other = -other
    end
    z = nil
    loop do
      z = z ? z * x : x if other[0] == 1
      return z if (other >>= 1).zero?
      x *= x
    end
  when Numeric
    v, d, v_inv = eigensystem
    v * Matrix.diagonal(*d.each(:diagonal).map{|e| e ** other}) * v_inv
  else
    Matrix.Raise ErrOperationNotDefined, "**", self.class, other.class
  end
end

# File ../ruby/lib/matrix.rb, line 835
def +(m)
  case m
  when Numeric
    Matrix.Raise ErrOperationNotDefined, "+", self.class, m.class
  when Vector
    m = Matrix.column_vector(m)
  when Matrix
  else
    return apply_through_coercion(m, __method__)
  end

  Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size

  rows = Array.new(row_size) {|i|
    Array.new(column_size) {|j|
      self[i, j] + m[i, j]
    }
  }
  new_matrix rows, column_size
end

# File ../ruby/lib/matrix.rb, line 862
def -(m)
  case m
  when Numeric
    Matrix.Raise ErrOperationNotDefined, "-", self.class, m.class
  when Vector
    m = Matrix.column_vector(m)
  when Matrix
  else
    return apply_through_coercion(m, __method__)
  end

  Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size

  rows = Array.new(row_size) {|i|
    Array.new(column_size) {|j|
      self[i, j] - m[i, j]
    }
  }
  new_matrix rows, column_size
end

# File ../ruby/lib/matrix.rb, line 889
def /(other)
  case other
  when Numeric
    rows = @rows.collect {|row|
      row.collect {|e| e / other }
    }
    return new_matrix rows, column_size
  when Matrix
    return self * other.inverse
  else
    return apply_through_coercion(other, __method__)
  end
end

# File ../ruby/lib/matrix.rb, line 764
def ==(other)
  return false unless Matrix === other &&
                      column_size == other.column_size # necessary for empty matrices
  rows == other.rows
end

# File ../ruby/lib/matrix.rb, line 313
def [](i, j)
  @rows.fetch(i){return nil}[j]
end

# File ../ruby/lib/matrix.rb, line 781
def clone
  new_matrix @rows.map(&:dup), column_size
end

# File ../ruby/lib/matrix.rb, line 1274
def coerce(other)
  case other
  when Numeric
    return Scalar.new(other), self
  else
    raise TypeError, "#{self.class} can't be coerced into #{other.class}"
  end
end

# File ../ruby/lib/matrix.rb, line 379
def collect(&block) # :yield: e
  return to_enum(:collect) unless block_given?
  rows = @rows.collect{|row| row.collect(&block)}
  new_matrix rows, column_size
end

# File ../ruby/lib/matrix.rb, line 356
def column(j) # :yield: e
  if block_given?
    return self if j >= column_size || j < -column_size
    row_size.times do |i|
      yield @rows[i][j]
    end
    self
  else
    return nil if j >= column_size || j < -column_size
    col = Array.new(row_size) {|i|
      @rows[i][j]
    }
    Vector.elements(col, false)
  end
end

# File ../ruby/lib/matrix.rb, line 1295
def column_vectors
  Array.new(column_size) {|i|
    column(i)
  }
end

# File ../ruby/lib/matrix.rb, line 1220
def conjugate
  collect(&:conjugate)
end

# File ../ruby/lib/matrix.rb, line 1006
def determinant
  Matrix.Raise ErrDimensionMismatch unless square?
  m = @rows
  case row_size
    # Up to 4x4, give result using Laplacian expansion by minors.
    # This will typically be faster, as well as giving good results
    # in case of Floats
  when 0
    +1
  when 1
    + m[0][0]
  when 2
    + m[0][0] * m[1][1] - m[0][1] * m[1][0]
  when 3
    m0, m1, m2 = m
    + m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1]        - m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0]        + m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0]
  when 4
    m0, m1, m2, m3 = m
    + m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2]        - m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1]        + m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1]        - m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2]        + m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0]        - m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0]        + m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1]        - m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0]        + m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0]        - m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1]        + m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0]        - m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0]
  else
    # For bigger matrices, use an efficient and general algorithm.
    # Currently, we use the Gauss-Bareiss algorithm
    determinant_bareiss
  end
end

# File ../ruby/lib/matrix.rb, line 1088
def determinant_e
  warn "#{caller(1)[0]}: warning: Matrix#determinant_e is deprecated; use #determinant"
  rank
end

# File ../ruby/lib/matrix.rb, line 597
def diagonal?
  Matrix.Raise ErrDimensionMismatch unless square?
  each(:off_diagonal).all?(&:zero?)
end

# File ../ruby/lib/matrix.rb, line 402
def each(which = :all) # :yield: e
  return to_enum :each, which unless block_given?
  last = column_size - 1
  case which
  when :all
    block = Proc.new
    @rows.each do |row|
      row.each(&block)
    end
  when :diagonal
    @rows.each_with_index do |row, row_index|
      yield row.fetch(row_index){return self}
    end
  when :off_diagonal
    @rows.each_with_index do |row, row_index|
      column_size.times do |col_index|
        yield row[col_index] unless row_index == col_index
      end
    end
  when :lower
    @rows.each_with_index do |row, row_index|
      0.upto([row_index, last].min) do |col_index|
        yield row[col_index]
      end
    end
  when :strict_lower
    @rows.each_with_index do |row, row_index|
      [row_index, column_size].min.times do |col_index|
        yield row[col_index]
      end
    end
  when :strict_upper
    @rows.each_with_index do |row, row_index|
      (row_index+1).upto(last) do |col_index|
        yield row[col_index]
      end
    end
  when :upper
    @rows.each_with_index do |row, row_index|
      row_index.upto(last) do |col_index|
        yield row[col_index]
      end
    end
  else
    Matrix.Raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
  end
  self
end

# File ../ruby/lib/matrix.rb, line 463
def each_with_index(which = :all) # :yield: e, row, column
  return to_enum :each_with_index, which unless block_given?
  last = column_size - 1
  case which
  when :all
    @rows.each_with_index do |row, row_index|
      row.each_with_index do |e, col_index|
        yield e, row_index, col_index
      end
    end
  when :diagonal
    @rows.each_with_index do |row, row_index|
      yield row.fetch(row_index){return self}, row_index, row_index
    end
  when :off_diagonal
    @rows.each_with_index do |row, row_index|
      column_size.times do |col_index|
        yield row[col_index], row_index, col_index unless row_index == col_index
      end
    end
  when :lower
    @rows.each_with_index do |row, row_index|
      0.upto([row_index, last].min) do |col_index|
        yield row[col_index], row_index, col_index
      end
    end
  when :strict_lower
    @rows.each_with_index do |row, row_index|
      [row_index, column_size].min.times do |col_index|
        yield row[col_index], row_index, col_index
      end
    end
  when :strict_upper
    @rows.each_with_index do |row, row_index|
      (row_index+1).upto(last) do |col_index|
        yield row[col_index], row_index, col_index
      end
    end
  when :upper
    @rows.each_with_index do |row, row_index|
      row_index.upto(last) do |col_index|
        yield row[col_index], row_index, col_index
      end
    end
  else
    Matrix.Raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
  end
  self
end

# File ../ruby/lib/matrix.rb, line 1187
def eigensystem
  EigenvalueDecomposition.new(self)
end

# File ../ruby/lib/matrix.rb, line 1308
def elements_to_f
  warn "#{caller(1)[0]}: warning: Matrix#elements_to_f is deprecated, use map(&:to_f)"
  map(&:to_f)
end

# File ../ruby/lib/matrix.rb, line 1313
def elements_to_i
  warn "#{caller(1)[0]}: warning: Matrix#elements_to_i is deprecated, use map(&:to_i)"
  map(&:to_i)
end

# File ../ruby/lib/matrix.rb, line 1318
def elements_to_r
  warn "#{caller(1)[0]}: warning: Matrix#elements_to_r is deprecated, use map(&:to_r)"
  map(&:to_r)
end

# File ../ruby/lib/matrix.rb, line 606
def empty?
  column_size == 0 || row_size == 0
end

# File ../ruby/lib/matrix.rb, line 770
def eql?(other)
  return false unless Matrix === other &&
                      column_size == other.column_size # necessary for empty matrices
  rows.eql? other.rows
end

# File ../ruby/lib/matrix.rb, line 788
def hash
  @rows.hash
end

# File ../ruby/lib/matrix.rb, line 614
def hermitian?
  Matrix.Raise ErrDimensionMismatch unless square?
  each_with_index(:strict_upper).all? do |e, row, col|
    e == rows[col][row].conj
  end
end

# File ../ruby/lib/matrix.rb, line 1234
def imaginary
  collect(&:imaginary)
end

# File ../ruby/lib/matrix.rb, line 526
def index(*args)
  raise ArgumentError, "wrong number of arguments(#{args.size} for 0-2)" if args.size > 2
  which = (args.size == 2 || SELECTORS.include?(args.last)) ? args.pop : :all
  return to_enum :find_index, which, *args unless block_given? || args.size == 1
  if args.size == 1
    value = args.first
    each_with_index(which) do |e, row_index, col_index|
      return row_index, col_index if e == value
    end
  else
    each_with_index(which) do |e, row_index, col_index|
      return row_index, col_index if yield e
    end
  end
  nil
end

# File ../ruby/lib/matrix.rb, line 1343
def inspect
  if empty?
    "Matrix.empty(#{row_size}, #{column_size})"
  else
    "Matrix#{@rows.inspect}"
  end
end

# File ../ruby/lib/matrix.rb, line 909
def inverse
  Matrix.Raise ErrDimensionMismatch unless square?
  Matrix.I(row_size).send(:inverse_from, self)
end

# File ../ruby/lib/matrix.rb, line 624
def lower_triangular?
  each(:strict_upper).all?(&:zero?)
end

# File ../ruby/lib/matrix.rb, line 1202
def lup
  LUPDecomposition.new(self)
end

# File ../ruby/lib/matrix.rb, line 556
def minor(*param)
  case param.size
  when 2
    row_range, col_range = param
    from_row = row_range.first
    from_row += row_size if from_row < 0
    to_row = row_range.end
    to_row += row_size if to_row < 0
    to_row += 1 unless row_range.exclude_end?
    size_row = to_row - from_row

    from_col = col_range.first
    from_col += column_size if from_col < 0
    to_col = col_range.end
    to_col += column_size if to_col < 0
    to_col += 1 unless col_range.exclude_end?
    size_col = to_col - from_col
  when 4
    from_row, size_row, from_col, size_col = param
    return nil if size_row < 0 || size_col < 0
    from_row += row_size if from_row < 0
    from_col += column_size if from_col < 0
  else
    Matrix.Raise ArgumentError, param.inspect
  end

  return nil if from_row > row_size || from_col > column_size || from_row < 0 || from_col < 0
  rows = @rows[from_row, size_row].collect{|row|
    row[from_col, size_col]
  }
  new_matrix rows, [column_size - from_col, size_col].min
end

# File ../ruby/lib/matrix.rb, line 632
def normal?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row_i, i|
    rows.each_with_index do |row_j, j|
      s = 0
      rows.each_with_index do |row_k, k|
        s += row_i[k] * row_j[k].conj - row_k[i].conj * row_k[j]
      end
      return false unless s == 0
    end
  end
  true
end

# File ../ruby/lib/matrix.rb, line 650
def orthogonal?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row, i|
    column_size.times do |j|
      s = 0
      row_size.times do |k|
        s += row[k] * rows[k][j]
      end
      return false unless s == (i == j ? 1 : 0)
    end
  end
  true
end

# File ../ruby/lib/matrix.rb, line 668
def permutation?
  Matrix.Raise ErrDimensionMismatch unless square?
  cols = Array.new(column_size)
  rows.each_with_index do |row, i|
    found = false
    row.each_with_index do |e, j|
      if e == 1
        return false if found || cols[j]
        found = cols[j] = true
      elsif e != 0
        return false
      end
    end
    return false unless found
  end
  true
end

# File ../ruby/lib/matrix.rb, line 1103
def rank
  # We currently use Bareiss' multistep integer-preserving gaussian elimination
  # (see comments on determinant)
  a = to_a
  last_column = column_size - 1
  last_row = row_size - 1
  pivot_row = 0
  previous_pivot = 1
  0.upto(last_column) do |k|
    switch_row = (pivot_row .. last_row).find {|row|
      a[row][k] != 0
    }
    if switch_row
      a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row
      pivot = a[pivot_row][k]
      (pivot_row+1).upto(last_row) do |i|
         ai = a[i]
         (k+1).upto(last_column) do |j|
           ai[j] =  (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot
         end
       end
      pivot_row += 1
      previous_pivot = pivot
    end
  end
  pivot_row
end

# File ../ruby/lib/matrix.rb, line 1134
def rank_e
  warn "#{caller(1)[0]}: warning: Matrix#rank_e is deprecated; use #rank"
  rank
end

# File ../ruby/lib/matrix.rb, line 1248
def real
  collect(&:real)
end

# File ../ruby/lib/matrix.rb, line 689
def real?
  all?(&:real?)
end

# File ../ruby/lib/matrix.rb, line 1258
def rect
  [real, imag]
end

# File ../ruby/lib/matrix.rb, line 696
def regular?
  not singular?
end

# File ../ruby/lib/matrix.rb, line 1142
def round(ndigits=0)
  map{|e| e.round(ndigits)}
end

# File ../ruby/lib/matrix.rb, line 342
def row(i, &block) # :yield: e
  if block_given?
    @rows.fetch(i){return self}.each(&block)
    self
  else
    Vector.elements(@rows.fetch(i){return nil})
  end
end

# File ../ruby/lib/matrix.rb, line 329
def row_size
  @rows.size
end

# File ../ruby/lib/matrix.rb, line 1286
def row_vectors
  Array.new(row_size) {|i|
    row(i)
  }
end

# File ../ruby/lib/matrix.rb, line 703
def singular?
  determinant == 0
end

# File ../ruby/lib/matrix.rb, line 710
def square?
  column_size == row_size
end

# File ../ruby/lib/matrix.rb, line 718
def symmetric?
  Matrix.Raise ErrDimensionMismatch unless square?
  each_with_index(:strict_upper).all? do |e, row, col|
    e == rows[col][row]
  end
end

# File ../ruby/lib/matrix.rb, line 1304
def to_a
  @rows.collect(&:dup)
end

# File ../ruby/lib/matrix.rb, line 1330
def to_s
  if empty?
    "Matrix.empty(#{row_size}, #{column_size})"
  else
    "Matrix[" + @rows.collect{|row|
      "[" + row.collect{|e| e.to_s}.join(", ") + "]"
    }.join(", ")+"]"
  end
end

# File ../ruby/lib/matrix.rb, line 1151
def trace
  Matrix.Raise ErrDimensionMismatch unless square?
  (0...column_size).inject(0) do |tr, i|
    tr + @rows[i][i]
  end
end

# File ../ruby/lib/matrix.rb, line 1169
def transpose
  return Matrix.empty(column_size, 0) if row_size.zero?
  new_matrix @rows.transpose, row_size
end

# File ../ruby/lib/matrix.rb, line 729
def unitary?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row, i|
    column_size.times do |j|
      s = 0
      row_size.times do |k|
        s += row[k].conj * rows[k][j]
      end
      return false unless s == (i == j ? 1 : 0)
    end
  end
  true
end

# File ../ruby/lib/matrix.rb, line 746
def upper_triangular?
  each(:strict_lower).all?(&:zero?)
end

# File ../ruby/lib/matrix.rb, line 753
def zero?
  all?(&:zero?)
end

# File ../ruby/lib/matrix.rb, line 319
def []=(i, j, v)
  @rows[i][j] = v
end

# File ../ruby/lib/matrix.rb, line 1057
def determinant_bareiss
  size = row_size
  last = size - 1
  a = to_a
  no_pivot = Proc.new{ return 0 }
  sign = +1
  pivot = 1
  size.times do |k|
    previous_pivot = pivot
    if (pivot = a[k][k]) == 0
      switch = (k+1 ... size).find(no_pivot) {|row|
        a[row][k] != 0
      }
      a[switch], a[k] = a[k], a[switch]
      pivot = a[k][k]
      sign = -sign
    end
    (k+1).upto(last) do |i|
      ai = a[i]
      (k+1).upto(last) do |j|
        ai[j] =  (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot
      end
    end
  end
  sign * pivot
end