Class: NMatrix

Inherits:
Object show all
Defined in:
lib/nmatrix/lapack.rb,
lib/nmatrix/nmatrix.rb,
lib/nmatrix/version.rb,
lib/nmatrix/shortcuts.rb

Overview

NMatrix

A linear algebra library for scientific computation in Ruby. NMatrix is part of SciRuby.

NMatrix was originally inspired by and derived from NArray, by Masahiro Tanaka: narray.rubyforge.org

Copyright Information

SciRuby is Copyright © 2010 - 2012, Ruby Science Foundation NMatrix is Copyright © 2012, Ruby Science Foundation

Please see LICENSE.txt for additional copyright notices.

Contributing

By contributing source code to SciRuby, you agree to be bound by our Contributor Agreement:

shortcuts.rb

These are shortcuts for NMatrix and NVector creation, contributed by Daniel Carrera ([email protected]) and Carlos Agarie ([email protected]).

Direct Known Subclasses

NVector

Defined Under Namespace

Modules: BLAS, IO, LAPACK

Constant Summary collapse

VERSION = 
'0.0.3'

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.bindgen(n) ⇒ Object 



224
225
226
 
# File 'lib/nmatrix/shortcuts.rb', line 224

def bindgen(n)
  NMatrix.seq(n, :byte)
end

.cindgen(n) ⇒ Object 



228
229
230
 
# File 'lib/nmatrix/shortcuts.rb', line 228

def cindgen(n)
  NMatrix.seq(n, :complex64)
end

.eye(*params) ⇒ Object Also known as: identity

identity() or eye()

Creates an identity matrix (square matrix rank 2) of the size
supplied as a parameter. Optional parameters include:

* A storage type as the first parameter (default is :dense).
* A dtype as the last parameter (default is :float64).

Examples:

eye(3) # =>   1.0   0.0   0.0
              0.0   1.0   0.0
              0.0   0.0   1.0

eye(3, :int32) # =>   1   0   0
                      0   1   0
                      0   0   1

eye(:yale, 2, :int32) # =>   1   0
                             0   1



103
104
105
106
107
108
109
110
111
112
113
114
115
116
 
# File 'lib/nmatrix/shortcuts.rb', line 103

def eye(*params)
  dtype = params.last.is_a?(Symbol) ? params.pop : :float64
  stype = params.first.is_a?(Symbol) ? params.shift : :dense

  dim = params.first
  
  # Fill the diagonal with 1's.
  m = NMatrix.zeros(stype, dim, dtype)
  (0 .. (dim - 1)).each do |i| 
    m[i, i] = 1
  end
  
  m
end

.findgen(n) ⇒ Object 



220
221
222
 
# File 'lib/nmatrix/shortcuts.rb', line 220

def findgen(n)
  NMatrix.seq(n, :float32)
end

.indgen(n) ⇒ Object

indgen() , findgen() , bindgen() , cindgen()

These IDL functions are similar to seq() but less flexible.
They produce one-dimensional vectors:

indgen    -- Integer vector   --   seq(n, :int32)
findgen   -- Float vector     --   seq(n, :float32)
bindgen   -- Byte vector      --   seq(n, :byte)
cindgen   -- Complex vector   --   seq(n, :complex64)



216
217
218
 
# File 'lib/nmatrix/shortcuts.rb', line 216

def indgen(n)
  NMatrix.seq(n, :int32)
end

.load_file(file_path) ⇒ Object 



198
199
200
 
# File 'lib/nmatrix/nmatrix.rb', line 198

def load_file(file_path)
	NMatrix::IO::Mat5Reader.new(File.open(file_path, 'rb')).to_ruby
end

.ones(*params) ⇒ Object

ones()

Creates a :dense matrix of ones with the dimensions supplied
as parameters. Optionaly, one can specify a dtype as the last
parameter (default is :float64).

Examples:

ones([1, 3]) # =>  1.0   1.0   1.0

ones([2, 3], :int32) # =>  1  1  1
                           1  1  1



74
75
76
77
78
79
 
# File 'lib/nmatrix/shortcuts.rb', line 74

def ones(*params)
  dtype = params.last.is_a?(Symbol) ? params.pop : :float64
  dim = params.first

  NMatrix.new(dim, 1, dtype)
end

.random(*params) ⇒ Object

random()

Creates a :dense NMatrix with random numbers between 0 and 1 generated
by Random::rand. The parameter is the dimension of the matrix.

Examples:

rand([2, 2]) # => 0.4859439730644226   0.1783195585012436
                  0.23193766176700592  0.4503345191478729



131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
 
# File 'lib/nmatrix/shortcuts.rb', line 131

def random(*params)
  dim = params.first
  rng = Random.new

  # Must provide the dimension as an Integer for a square matrix or as an
  # array, e.g. [2, 4, 7].
  unless dim.is_a?(Integer) || dim.is_a?(Array)
    raise ArgumentError, "random() accepts only integers or arrays as \
dimension."
  end

  random_values = []
  
  # Construct the values of the final matrix based on the dimension.
  if dim.is_a?(Integer)
    (dim * dim - 1).times { |i| random_values << rng.rand }
  else
    # Dimensions given by an array. Get the product of the array elements
    # and generate this number of random values.
    dim.reduce(1, :*).times { |i| random_values << rng.rand }
  end

  NMatrix.new(:dense, dim, random_values, :float64)
end

.seq(*params) ⇒ Object

seq()

Creates a :dense NMatrix with a sequence of integers starting at
zero until the matrix is filled. The parameters to the method
are the dimensions of the matrix. Optionaly, one can specify a
dtype as the last parameter (default is :float64).

Examples:

seq(2) # =>   0   1
              2   3

seq([3, 3], :float32) # =>  0.0  1.0  2.0
                            3.0  4.0  5.0
                            6.0  7.0  8.0



173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
 
# File 'lib/nmatrix/shortcuts.rb', line 173

def seq(*params)
  dtype = params.last.is_a?(Symbol) ? params.pop : nil
  dim = params.first
  
  # Must provide the dimension as an Integer for a square matrix or as an
  # 2 element array, e.g. [2,4].
  unless dim.is_a?(Integer) || (dim.is_a?(Array) && dim.size < 3)
    raise ArgumentError, "seq() accepts only integers or 2-element arrays \
as dimension."
  end
  
  # Construct the values of the final matrix based on the dimension.
  if dim.is_a?(Integer)
    values = (0 .. (dim * dim - 1)).to_a
  else
    # Dimensions given by a 2 element array.
    values = (0 .. (dim.first * dim.last - 1)).to_a
  end
  
  # It'll produce :int32, except if a dtype is provided.
  NMatrix.new(:dense, dim, values, dtype)
end

.zeros(*params) ⇒ Object Also known as: zeroes

zeros() or zeroes()

Creates a new matrix of zeros with the dimensions supplied as
parameters. Optional parameters include:

* A storage type as the first parameter (default is :dense).
* A dtype as the last parameter (default is :float64).

Examples:

zeros(2) # =>  0.0   0.0   
               0.0   0.0

zeros([2, 3], :int32) # =>  0  0  0
                         0  0  0

zeros(:list, [1, 5], :int32) # =>  0  0  0  0  0



50
51
52
53
54
55
56
 
# File 'lib/nmatrix/shortcuts.rb', line 50

def zeros(*params)
  dtype = params.last.is_a?(Symbol) ? params.pop : :float64
  stype = params.first.is_a?(Symbol) ? params.shift : :dense
  dim = params.first

  NMatrix.new(stype, dim, 0, dtype)
end

Instance Method Details

#__yale_ary__to_s(sym) ⇒ Object 



191
192
193
194
195
 
# File 'lib/nmatrix/nmatrix.rb', line 191

def __yale_ary__to_s(sym)
	ary = self.send("__yale_#{sym.to_s}__".to_sym)
	
	'[' + ary.collect { |a| a ? a : 'nil'}.join(',') + ']'
end

#cols ⇒ Object 



94
95
96
 
# File 'lib/nmatrix/nmatrix.rb', line 94

def cols
  shape[1]
end

#column(column_number, get_by = :copy) ⇒ Object

column()

Returns the column specified. The second parameter defaults to
:copy, which returns a copy of the selected column, but it can be
specified as :reference, which will return a reference to it.

Examples:

m = NMatrix.new(2, [1, 4, 9, 14], :int32) # =>  1   4
                                                9  14

m.column(1) # =>   4
                  14



256
257
258
259
260
261
262
263
264
265
266
 
# File 'lib/nmatrix/shortcuts.rb', line 256

def column(column_number, get_by = :copy)
  unless [:copy, :reference].include?(get_by)
    raise ArgumentError, "column() 2nd parameter must be :copy or :reference"
  end
  
  if get_by == :copy
    self.slice(0 ... self.shape[0], column_number)
  else # by reference
    self[0 ... self.shape[0], column_number]
  end
end

#complex_conjugate(new_stype = self.stype) ⇒ Object

Get the complex conjugate of this matrix. See also complex_conjugate! for an in-place operation (provided the dtype is already :complex64 or :complex128).

Does not work on list matrices, but you can optionally pass in the type you want to cast to if you’re dealing with a list matrix.



164
165
166
 
# File 'lib/nmatrix/nmatrix.rb', line 164

def complex_conjugate(new_stype = self.stype)
	self.cast(new_stype, NMatrix::upcast(dtype, :complex64)).complex_conjugate!
end

#conjugate_transpose ⇒ Object

Calculate the conjugate transpose of a matrix. If your dtype is already complex, this should only require one copy (for the transpose).



170
171
172
 
# File 'lib/nmatrix/nmatrix.rb', line 170

def conjugate_transpose
	self.transpose.complex_conjugate!
end

#det ⇒ Object

Calculate the determinant by way of LU decomposition. This is accomplished using clapack_getrf, and then by summing the diagonal elements. There is a risk of underflow/overflow.

There are probably also more efficient ways to calculate the determinant. This method requires making a copy of the matrix, since clapack_getrf modifies its input.

For smaller matrices, you may be able to use det_exact.

This function is guaranteed to return the same type of data in the matrix upon which it is called. In other words, if you call it on a rational matrix, you’ll get a rational number back.

Integer matrices are converted to rational matrices for the purposes of performing the calculation, as xGETRF can’t work on integer matrices.

Raises:

  • (NotImplementedError) 


138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
 
# File 'lib/nmatrix/nmatrix.rb', line 138

def det
  raise(NotImplementedError, "determinant can be calculated only for 2D matrices") unless self.dim == 2

  # Cast to a dtype for which getrf is implemented
  new_dtype = [:byte,:int8,:int16,:int32,:int64].include?(self.dtype) ? :rational128 : self.dtype
  copy = self.cast(:dense, new_dtype)

  # Need to know the number of permutations. We'll add up the diagonals of the factorized matrix.
  pivot = copy.getrf!

  prod = pivot.size % 2 == 1 ? -1 : 1 # odd permutations => negative
  [shape[0],shape[1]].min.times do |i|
    prod *= copy[i,i]
  end

  # Convert back to an integer if necessary
  new_dtype != self.dtype ? prod.to_i : prod
end

#getrf! ⇒ Object

Calls clapack_getrf and returns the pivot array (dense only).

Raises:

  • (StorageTypeError) 


119
120
121
122
 
# File 'lib/nmatrix/nmatrix.rb', line 119

def getrf!
  raise(StorageTypeError, "ATLAS functions only work on dense matrices") unless self.stype == :dense
  NMatrix::LAPACK::clapack_getrf(:row, self.shape[0], self.shape[1], self, self.shape[0])
end

#hermitian? ⇒ Boolean

Returns:

  • (Boolean) 


174
175
176
177
178
179
180
181
182
183
 
# File 'lib/nmatrix/nmatrix.rb', line 174

def hermitian?
	return false if self.dim != 2 or self.shape[0] != self.shape[1]
	
	if [:complex64, :complex128].include?(self.dtype)
		# TODO: Write much faster Hermitian test in C
		self.eql?(self.conjugate_transpose)
	else
		symmetric?
	end
end

#inspect ⇒ Object 



185
186
187
188
189
 
# File 'lib/nmatrix/nmatrix.rb', line 185

def inspect
	original_inspect = super()
	original_inspect = original_inspect[0...original_inspect.size-1]
	original_inspect + inspect_helper.join(" ") + ">"
end

#invert ⇒ Object Also known as: inverse

Make a copy of the matrix, then invert it (requires LAPACK). Returns a dense matrix.



112
113
114
 
# File 'lib/nmatrix/nmatrix.rb', line 112

def invert
  self.cast(:dense, self.dtype).invert!
end

#invert! ⇒ Object

Use LAPACK to calculate the inverse of the matrix (in-place). Only works on dense matrices.

Note: If you don’t have LAPACK, e.g., on a Mac, this may not work yet.



101
102
103
104
105
106
107
108
109
 
# File 'lib/nmatrix/nmatrix.rb', line 101

def invert!
  # Get the pivot array; factor the matrix
  pivot = self.getrf!

  # Now calculate the inverse using the pivot array
  NMatrix::LAPACK::clapack_getri(:row, self.shape[0], self, self.shape[0], pivot)

  self
end

#pretty_print(q = nil) ⇒ Object Also known as: pp

TODO: Make this actually pretty.



59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
 
# File 'lib/nmatrix/nmatrix.rb', line 59

def pretty_print(q = nil)
	if dim != 2 || (dim == 2 && shape[1] > 10) # FIXME: Come up with a better way of restricting the display
     inspect
   else

     arr = (0...shape[0]).map do |i|
       ary = []
       (0...shape[1]).each do |j|
         o = begin
           self[i, j]
         rescue ArgumentError
           nil
         end
         ary << (o.nil? ? 'nil' : o)
       end
       ary.inspect
     end

     if q.nil?
       puts arr.join("\n")
     else
       q.group(1, "", "\n") do
         q.seplist(arr, lambda { q.text "  " }, :each)  { |v| q.text v.to_s }
       end
     end

   end
end

#rows ⇒ Object

These shortcuts use #shape to return the number of rows and columns.



90
91
92
 
# File 'lib/nmatrix/nmatrix.rb', line 90

def rows
  shape[0]
end