class Matrix

# File lib/extendmatrix.rb, line 320
def build(n,m,val=0)
  f = (block_given?)? lambda {|i,j| yield(i, j)} : lambda {|i,j| val}
  Matrix.rows((0...n).collect {|i| (0...m).collect {|j| f.call(i,j)}})
end

# File lib/extendmatrix.rb, line 328
def diag(*args)
  dsize = 0
  sizes = args.collect{|e| x = (e.is_a?(Matrix)) ? e.row_size : 1; dsize += x; x}
  m = Matrix.zero(dsize)
  count = 0

  sizes.size.times{|i|
    range = count..(count+sizes[i]-1)
    m[range, range] = args[i]
    count += sizes[i]
  }
  m
end

# File lib/extendmatrix.rb, line 359
def diag_in_delta?(m1, m0, eps = 1.0e-10)
  n = m1.row_size
  return false if n != m0.row_size or m1.column_size != m0.column_size
  n.times{|i|
    return false if (m1[i,i]-m0[i,i]).abs > eps
  }
  true
end

# File lib/extendmatrix.rb, line 345
def equal_in_delta?(m0, m1, delta = 1.0e-10)
  delta = delta.abs
  m0.row_size.times {|i|
    m0.column_size.times {|j|
      x=m0[i,j]; y=m1[i,j]
      return false if (x < y - delta or x > y + delta)
    }
  }
  true
end

# File lib/extendmatrix.rb, line 226
def self.robust_I(n)
  n>0 ? self.I(n) : self.[]
end

# File lib/extendmatrix.rb, line 241
def [](i, j)
  case i
  when Range
    case j
    when Range
      Matrix[*i.collect{|l| self.row(l)[j].to_a}]
    else
      column(j)[i]
    end
  else
    case j
    when Range
      row(i)[j]
    else
      ids(i, j)
    end
  end
end

# File lib/extendmatrix.rb, line 270
def []=(i, j, v)
  case i
  when Range
    if i.entries.size == 1
      self[i.begin, j] = (v.is_a?(Matrix) ? v.row(0) : v)
    else
      case j
      when Range
        if j.entries.size == 1
          self[i, j.begin] = (v.is_a?(Matrix) ? v.column(0) : v)
        else
          i.each{|l| self.row= l, v.row(l - i.begin), j}
        end
      else
        self.column= j, v, i
      end
    end
  else
    case j
    when Range
      if j.entries.size == 1
        self[i, j.begin] = (v.is_a?(Vector) ? v[0] : v)
      else
        self.row= i, v, j
      end
    else
      @rows[i][j] = (v.is_a?(Vector) ? v[0] : v)

    end
  end
end

# File lib/extendmatrix.rb, line 429
def cols_len
  (0...column_size).collect {|j| max_len_column(j)}
end

# File lib/extendmatrix.rb, line 304
def dup
  (self.class).rows(self.rows.dup)
end

# File lib/extendmatrix.rb, line 463
def each
  @rows.each {|x| x.each {|e| yield(e)}}
  nil
end

# File lib/extendmatrix.rb, line 308
def initialize_copy(orig)
  init_rows(orig.rows, true)
  self.wrap=(orig.wrap)
end

# File lib/extendmatrix.rb, line 422
def max_len_column(j)
  column_collect(j) {|x| x.to_s.length}.max
end

# File lib/extendmatrix.rb, line 372
def quo(v)
  map {|e| e.quo(v)}
end

# File lib/extendmatrix.rb, line 385
def set(m)
  0.upto(m.row_size - 1) do |i|
    0.upto(m.column_size - 1) do |j|
      self[i, j] = m[i, j]
    end
  end
  self.wrap = m.wrap
end

# File lib/extendmatrix.rb, line 437
def to_s(mode = :pretty, len_col = 3)
  return "Matrix[]" if empty?
  if mode == :pretty
    clen = cols_len
    to_a.collect {|r|
      i=0
      r.map {|x|
        l=clen[i]
        i+=1
        format("%#{l}s ", x.to_s)
      } << "\n"
    }.join("")
  else
    i = 0; s = ""; cs = column_size
    each do |e|
      i = (i + 1) % cs
      s += format("%#{len_col}s ", e.to_s)
      s += "\n" if i == 0
    end
    s
  end
end

# File lib/extendmatrix.rb, line 409
def wrap=(mode = :torus)
  case mode
  when :torus then eval(wraplate("i %= row_size; j %= column_size"))
  when :h_cylinder then eval(wraplate("i %= row_size"))
  when :v_cylinder then eval(wraplate("j %= column_size"))
  when :nil then eval(wraplate)
  end
  @wrap = mode
end

# File lib/extendmatrix.rb, line 394
def wraplate(ijwrap = "")
  "class << self
    def [](i, j)
    #{ijwrap}; @rows[i][j]
    end

    def []=(i, j, v)
    #{ijwrap}; @rows[i][j] = v
    end
  end"
end

# File lib/extendmatrix.rb, line 744
def L
  LU.factorization(self)[0]
end

# File lib/extendmatrix.rb, line 737
def U
  LU.factorization(self)[1]
end

# File lib/extendmatrix.rb, line 825
def bidiagonal
  ub, vb = Householder.bidiag(self)
  ub.t * self * vb
end

# File lib/extendmatrix.rb, line 1074
def cJacobi(tol = 1.0e-10)
  a = self.clone
  n = row_size
  v = Matrix.I(n)
  eps = tol * a.normF
  while Jacobi.off(a) > eps
    p, q = Jacobi.max(a)
    c, s = Jacobi.sym_schur2(a, p, q)
    #print "\np:#{p} q:#{q} c:#{c} s:#{s}\n"
    j = Jacobi.J(p, q, c, s, n)
    a = j.t * a * j
    v = v * j
  end
  return a, v
end

# File lib/extendmatrix.rb, line 1094
def cJacobiA(tol = 1.0e-10)
  cJacobi(tol)[0]
end

# File lib/extendmatrix.rb, line 1111
def cJacobiV(tol = 1.0e-10)
  cJacobi(tol)[1]
end

# File lib/extendmatrix.rb, line 591
def column!(j)
  if block_given?
    (0...row_size).collect { |i| @rows[i][j] = yield @rows[i][j] }
  else
    column(j)
  end
end

# File lib/extendmatrix.rb, line 665
def column2matrix(c)
  if c.is_a?(Range)
    a=c.map {|v| column(v).to_a}.find_all {|v| v.size>0}
  else
    a = column(c).to_a
  end
  if c.is_a?(Range)
    return Matrix.columns(a)
  else
    return Matrix[*a.collect{|x| [x]}]
  end
end

# File lib/extendmatrix.rb, line 608
def column=(args)
  m = row_size
  c, v, r = MMatrix.id_vect_range(args, m)
  (m..r.begin - 1).each{|i| self[i, c] = 0}
  [v.size, r.entries.size].min.times{|i| self[i + r.begin, c] = v[i]}
  ((v.size + r.begin)..r.entries.last).each {|i| self[i, c] = 0}
end

# File lib/extendmatrix.rb, line 582
def column_collect(j, &block)
  f = MMatrix.default_block(block)
  (0...row_size).collect {|r| f.call(self[r, j])}
end

# File lib/extendmatrix.rb, line 685
def column_sum
  (0...column_size).collect {|i|
    column(i).sum
  }
end

# File lib/extendmatrix.rb, line 1102
def eigenvaluesJacobi
  a = cJacobiA
  Vector[*(0...row_size).collect{|i| a[i, i]}]
end

# File lib/extendmatrix.rb, line 641
def empty?
  @rows.empty? if @rows
end

# File lib/extendmatrix.rb, line 941
def givensQ
  Givens.QR(self)[1]
end

# File lib/extendmatrix.rb, line 932
def givensR
  Givens.QR(self)[0]
end

# File lib/extendmatrix.rb, line 856
def gram_schmidt
  a = clone
  n = column_size
  m = row_size
  q = Matrix.build(m, n){0}
  r = Matrix.zero(n)
  for k in 0...n
    r[k,k] = a[0...m, k].norm
    q[0...m, k] = a[0...m, k] / r[k, k]
    for j in (k+1)...n
      r[k, j] = q[0...m, k].t * a[0...m, j]
      a[0...m, j] -= q[0...m, k] * r[k, j]
    end
  end
  return q, r
end

# File lib/extendmatrix.rb, line 877
def gram_schmidtQ
  gram_schmidt[0]
end

# File lib/extendmatrix.rb, line 884
def gram_schmidtR
  gram_schmidt[1]
end

# File lib/extendmatrix.rb, line 967
def hessenbergQ
  Hessenberg.QR(self)[0]
end

# File lib/extendmatrix.rb, line 974
def hessenbergR
  Hessenberg.QR(self)[1]
end

# File lib/extendmatrix.rb, line 981
def hessenberg_form_H
  Householder.toHessenberg(self)[0]
end

# File lib/extendmatrix.rb, line 834
def houseQ
  h = Householder.QR(self)
  q = h[0]
  (1...h.size).each{|i| q *= h[i]}
  q
end

# File lib/extendmatrix.rb, line 845
def houseR
  h = Householder.QR(self)
  r = self.clone
  h.size.times{|i| r = h[i] * r}
  r
end

# File lib/extendmatrix.rb, line 629
def norm(p = 2)
  Vector::Norm.sqnorm(self, p) ** (1.quo(p))
end

# File lib/extendmatrix.rb, line 633
def norm_frobenius
  norm
end

# File lib/extendmatrix.rb, line 989
def realSchur(eps = 1.0e-10, steps = 100)
  h = self.hessenberg_form_H
  h1 = Matrix[]
  i = 0
  loop do
    h1 = h.hessenbergR * h.hessenbergQ
    break if Matrix.diag_in_delta?(h1, h, eps) or steps <= 0
    h = h1.clone
    steps -= 1
    i += 1
  end
  h1
end

# File lib/extendmatrix.rb, line 570
def row!(i)
  if block_given?
    @rows[i].collect! {|e| yield e }
  else
    Vector.elements(@rows[i], false)
  end
end

# File lib/extendmatrix.rb, line 649
def row2matrix(r)
  if r.is_a? Range
    a=r.map {|v| v<row_size ? row(v).to_a : nil }.find_all {|v| !v.nil?}
  else
    a = row(r).to_a
  end
  if r.is_a?(Range) and r.entries.size > 1
    return Matrix[*a]
  else
    return Matrix[a]
  end
end

# File lib/extendmatrix.rb, line 624
def row=(args)
  i, val, range = MMatrix.id_vect_range(args, column_size)
  row!(i)[range] = val
end

# File lib/extendmatrix.rb, line 562
def row_collect(i, &block)
  f = MMatrix.default_block(block)
  @rows[i].collect {|e| f.call(e)}
end

# File lib/extendmatrix.rb, line 679
def row_sum
  (0...row_size).collect {|i|
    row(i).sum
  }
end

# File lib/extendmatrix.rb, line 691
def total_sum
  row_sum.inject(&:+)
end

# File lib/extendmatrix.rb, line 526
def diagonal
  min = (row_size<column_size) ? row_size : column_size
  min.times.map {|i| self[i,i]}
end

# File lib/extendmatrix.rb, line 479
def e_mult(other)
  elementwise_operation(:*,other)
end

# File lib/extendmatrix.rb, line 484
def e_quo(other)
  elementwise_operation(:quo,other)
end

# File lib/extendmatrix.rb, line 505
def eigen
  ep=eigenpairs
  { :eigenvalues=>ep.map {|a| a[0]} ,
    :eigenvectors=>Matrix.columns(ep.map{|a| a[1]})
  }
end

# File lib/extendmatrix.rb, line 493
def eigenpairs
  eigval, eigvec= eigenvaluesJacobi, cJacobiV
  eigenpairs=eigval.size.times.map {|i|
    [eigval[i], eigvec.column(i)]
  }
  eigenpairs=eigenpairs.sort{|a,b| a[0]<=>b[0]}.reverse
end

# File lib/extendmatrix.rb, line 472
def elementwise_operation(op,other)
  Matrix.build(row_size,column_size) do |row, column|
    self[row,column].send(op,other[row,column])
  end
end

# File lib/extendmatrix.rb, line 512
def ln
  map {|v| Math::log(v)}
end

# File lib/extendmatrix.rb, line 489
def mssq
  @rows.inject(0){|ac,row| ac+(row.inject(0) {|acr,i| acr+(i**2)})}
end

# File lib/extendmatrix.rb, line 516
def sqrt
  map {|v| Math::sqrt(v)}
end

# File lib/extendmatrix.rb, line 520
def sscp
  transpose*self
end