Matrix

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

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

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

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

 
               # File 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 matrix.rb, line 287
def Matrix.empty(row_count = 0, column_count = 0)
  raise ArgumentError, "One size must be 0" if column_count != 0 && row_count != 0
  raise ArgumentError, "Negative size" if column_count < 0 || row_count < 0

  new([[]]*row_count, column_count)
end
            

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

 
               # File matrix.rb, line 297
def initialize(rows, column_count = rows[0].size)
  # No checking is done at this point. rows must be an Array of Arrays.
  # column_count 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_count = column_count
end
            

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

 
               # File 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|
    raise ErrDimensionMismatch, "row size differs (#{row.size} should be #{size})" unless row.size == size
  end
  new rows, size
end
            

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

 
               # File matrix.rb, line 245
def Matrix.zero(row_count, column_count = row_count)
  rows = Array.new(row_count){Array.new(column_count, 0)}
  new rows, column_count
end
            

 
               # File matrix.rb, line 805
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_count
  when Vector
    m = self.class.column_vector(m)
    r = self * m
    return r.column(0)
  when Matrix
    Matrix.Raise ErrDimensionMismatch if column_count != m.row_count

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

 
               # File matrix.rb, line 972
def ** (other)
  case other
  when Integer
    x = self
    if other <= 0
      x = self.inverse
      return self.class.identity(self.column_count) 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 * self.class.diagonal(*d.each(:diagonal).map{|e| e ** other}) * v_inv
  else
    Matrix.Raise ErrOperationNotDefined, "**", self.class, other.class
  end
end
            

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

  Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count and column_count == m.column_count

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

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

  Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count and column_count == m.column_count

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

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

 
               # File matrix.rb, line 767
def ==(other)
  return false unless Matrix === other &&
                      column_count == other.column_count # necessary for empty matrices
  rows == other.rows
end
            

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

 
               # File matrix.rb, line 784
def clone
  new_matrix @rows.map(&:dup), column_count
end
            

 
               # File matrix.rb, line 1277
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 matrix.rb, line 381
def collect(&block) # :yield: e
  return to_enum(:collect) unless block_given?
  rows = @rows.collect{|row| row.collect(&block)}
  new_matrix rows, column_count
end
            

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

 
               # File matrix.rb, line 1298
def column_vectors
  Array.new(column_count) {|i|
    column(i)
  }
end
            

 
               # File matrix.rb, line 1223
def conjugate
  collect(&:conjugate)
end
            

 
               # File matrix.rb, line 1009
def determinant
  Matrix.Raise ErrDimensionMismatch unless square?
  m = @rows
  case row_count
    # 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 matrix.rb, line 1091
def determinant_e
  warn "#{caller(1)[0]}: warning: Matrix#determinant_e is deprecated; use #determinant"
  determinant
end
            

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

 
               # File matrix.rb, line 404
def each(which = :all) # :yield: e
  return to_enum :each, which unless block_given?
  last = column_count - 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_count.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_count].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
    raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
  end
  self
end
            

 
               # File matrix.rb, line 465
def each_with_index(which = :all) # :yield: e, row, column
  return to_enum :each_with_index, which unless block_given?
  last = column_count - 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_count.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_count].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
    raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
  end
  self
end
            

 
               # File matrix.rb, line 1190
def eigensystem
  EigenvalueDecomposition.new(self)
end
            

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

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

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

 
               # File matrix.rb, line 608
def empty?
  column_count == 0 || row_count == 0
end
            

 
               # File matrix.rb, line 773
def eql?(other)
  return false unless Matrix === other &&
                      column_count == other.column_count # necessary for empty matrices
  rows.eql? other.rows
end
            

 
               # File matrix.rb, line 791
def hash
  @rows.hash
end
            

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

 
               # File matrix.rb, line 1237
def imaginary
  collect(&:imaginary)
end
            

 
               # File matrix.rb, line 528
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 matrix.rb, line 1346
def inspect
  if empty?
    "#{self.class}.empty(#{row_count}, #{column_count})"
  else
    "#{self.class}#{@rows.inspect}"
  end
end
            

 
               # File matrix.rb, line 912
def inverse
  Matrix.Raise ErrDimensionMismatch unless square?
  self.class.I(row_count).send(:inverse_from, self)
end
            

 
               # File matrix.rb, line 626
def lower_triangular?
  each(:strict_upper).all?(&:zero?)
end
            

 
               # File matrix.rb, line 1205
def lup
  LUPDecomposition.new(self)
end
            

 
               # File matrix.rb, line 558
def minor(*param)
  case param.size
  when 2
    row_range, col_range = param
    from_row = row_range.first
    from_row += row_count if from_row < 0
    to_row = row_range.end
    to_row += row_count 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_count if from_col < 0
    to_col = col_range.end
    to_col += column_count 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_count if from_row < 0
    from_col += column_count if from_col < 0
  else
    raise ArgumentError, param.inspect
  end

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

 
               # File matrix.rb, line 634
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 matrix.rb, line 652
def orthogonal?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row, i|
    column_count.times do |j|
      s = 0
      row_count.times do |k|
        s += row[k] * rows[k][j]
      end
      return false unless s == (i == j ? 1 : 0)
    end
  end
  true
end
            

 
               # File matrix.rb, line 670
def permutation?
  Matrix.Raise ErrDimensionMismatch unless square?
  cols = Array.new(column_count)
  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 matrix.rb, line 1106
def rank
  # We currently use Bareiss' multistep integer-preserving gaussian elimination
  # (see comments on determinant)
  a = to_a
  last_column = column_count - 1
  last_row = row_count - 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 matrix.rb, line 1137
def rank_e
  warn "#{caller(1)[0]}: warning: Matrix#rank_e is deprecated; use #rank"
  rank
end
            

 
               # File matrix.rb, line 1251
def real
  collect(&:real)
end
            

 
               # File matrix.rb, line 691
def real?
  all?(&:real?)
end
            

 
               # File matrix.rb, line 1261
def rect
  [real, imag]
end
            

 
               # File matrix.rb, line 698
def regular?
  not singular?
end
            

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

 
               # File matrix.rb, line 344
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 matrix.rb, line 329
def row_count
  @rows.size
end
            

 
               # File matrix.rb, line 1289
def row_vectors
  Array.new(row_count) {|i|
    row(i)
  }
end
            

 
               # File matrix.rb, line 705
def singular?
  determinant == 0
end
            

 
               # File matrix.rb, line 712
def square?
  column_count == row_count
end
            

 
               # File matrix.rb, line 720
def symmetric?
  Matrix.Raise ErrDimensionMismatch unless square?
  each_with_index(:strict_upper) do |e, row, col|
    return false if e != rows[col][row]
  end
  true
end
            

 
               # File matrix.rb, line 1307
def to_a
  @rows.collect(&:dup)
end
            

 
               # File matrix.rb, line 1333
def to_s
  if empty?
    "#{self.class}.empty(#{row_count}, #{column_count})"
  else
    "#{self.class}[" + @rows.collect{|row|
      "[" + row.collect{|e| e.to_s}.join(", ") + "]"
    }.join(", ")+"]"
  end
end
            

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

 
               # File matrix.rb, line 1172
def transpose
  return self.class.empty(column_count, 0) if row_count.zero?
  new_matrix @rows.transpose, row_count
end
            

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

 
               # File matrix.rb, line 749
def upper_triangular?
  each(:strict_lower).all?(&:zero?)
end
            

 
               # File matrix.rb, line 756
def zero?
  all?(&:zero?)
end