#--
# = 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: http://narray.rubyforge.org
#
# == Copyright Information
#
# SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
# NMatrix is Copyright (c) 2012 - 2014, John Woods and the 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:
#
# * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement
#
# == shortcuts.rb
#
# These are shortcuts for NMatrix and NVector creation, contributed by Daniel
# Carrera (dcarrera@hush.com) and Carlos Agarie (carlos.agarie@gmail.com).
#
# TODO Make all the shortcuts available through modules, allowing someone
# to include them to make "MATLAB-like" scripts.
#
# There are some questions to be answered before this can be done, tho.
#++
class NMatrix
# Methods for generating magic matrix.
module MagicHelpers
class << self
def odd_magic(nm, shape)
row = shape - 1
col = shape / 2
nm[row,col] = 1
(2..shape * shape).each do |index|
if nm[(row + 1) % shape,(col + 1) % shape] == 0
row = (row + 1) % shape
col = (col + 1) % shape
else
row = (row - 1 + shape) % shape
end
nm[row,col] = index
end
end
def doubly_even_magic(nm, shape)
mini_square_num = shape / 4
count = 1
inv_count = shape * shape
shape.times do |row|
shape.times do |col|
if col >= mini_square_num and col < shape - mini_square_num
if row >= mini_square_num and row < shape - mini_square_num
nm[row,col] = count
else
nm[row,col] = inv_count
end
elsif row < mini_square_num or row >= shape - mini_square_num
nm[row,col] = count
else
nm[row,col] = inv_count
end
count += 1
inv_count -= 1
end
end
end
def singly_even_magic(nm, shape)
half_shape = shape / 2
complementary_pair = (shape - 2) / 4
swap_col = NMatrix.new([shape])
index = 0
mini_magic = NMatrix.new([half_shape,half_shape], 0, dtype: nm.dtype)
odd_magic mini_magic, half_shape
half_shape.times do |row|
half_shape.times do |col|
nm[row,col] = mini_magic[row,col]
nm[row + half_shape,col + half_shape] = mini_magic[row,col] + half_shape * half_shape
nm[row,col + half_shape] = mini_magic[row,col] + 2 * half_shape * half_shape
nm[row + half_shape,col] = mini_magic[row,col] + 3 * half_shape * half_shape
end
end
(1..complementary_pair).each do |complementary_entry|
swap_col[index] = complementary_entry
index += 1
end
(shape - complementary_pair + 2..shape).each do |center|
swap_col[index] = center
index += 1
end
(1..half_shape).each do |row|
(1..index).each do |col|
temp = nm[row - 1,swap_col[col - 1] - 1]
nm[row - 1,swap_col[col - 1] - 1] = nm[row + half_shape - 1,swap_col[col - 1] - 1]
nm[row + half_shape - 1,swap_col[col - 1] - 1] = temp
end
end
temp = nm[complementary_pair,0]
nm[complementary_pair,0] = nm[complementary_pair + half_shape,0]
nm[complementary_pair + half_shape,0] = temp
temp = nm[complementary_pair + half_shape,complementary_pair]
nm[complementary_pair + half_shape,complementary_pair] = nm[complementary_pair,complementary_pair]
nm[complementary_pair,complementary_pair] = temp
end
end
end
# call-seq:
# m.dense? -> true or false
#
# Determine if +m+ is a dense matrix.
def dense?; return stype == :dense; end
# call-seq:
# m.yale? -> true or false
#
# Determine if +m+ is a Yale matrix.
def yale?; return stype == :yale; end
# call-seq:
# m.list? -> true or false
#
# Determine if +m+ is a list-of-lists matrix.
def list?; return stype == :list; end
class << self
# call-seq:
# NMatrix[Numeric, ..., Numeric, dtype: Symbol] -> NMatrix
# NMatrix[Array, dtype: Symbol] -> NMatrix
#
# The default value for +dtype+ is guessed from the first parameter. For example:
# NMatrix[1.0, 2.0].dtype # => :float64
#
# But this is just a *guess*. If the other values can't be converted to
# this dtype, a +TypeError+ will be raised.
#
# You can use the +N+ constant in this way:
# N = NMatrix
# N[1, 2, 3]
#
# NMatrix needs to have a succinct way to create a matrix by specifying the
# components directly. This is very useful for using it as an advanced
# calculator, it is useful for learning how to use, for testing language
# features and for developing algorithms.
#
# The NMatrix::[] method provides a way to create a matrix in a way that is compact and
# natural. The components are specified using Ruby array syntax. Optionally,
# one can specify a dtype as the last parameter (default is :float64).
#
# Examples:
#
# a = N[ 1,2,3,4 ] => 1 2 3 4
#
# a = N[ 1,2,3,4, :int32 ] => 1 2 3 4
#
# a = N[ [1,2,3], [3,4,5] ] => 1.0 2.0 3.0
# 3.0 4.0 5.0
#
# a = N[ 3,6,9 ].transpose => 3
# 6
# 9
#
# SYNTAX COMPARISON:
#
# MATLAB: a = [ [1 2 3] ; [4 5 6] ] or [ 1 2 3 ; 4 5 6 ]
# IDL: a = [ [1,2,3] , [4,5,6] ]
# NumPy: a = array( [1,2,3], [4,5,6] )
#
# SciRuby: a = NMatrix[ [1,2,3], [4,5,6] ]
# Ruby array: a = [ [1,2,3], [4,5,6] ]
def [](*params)
options = params.last.is_a?(Hash) ? params.pop : {}
# First find the dimensions of the array.
i = 0
shape = []
row = params
while row.is_a?(Array)
shape[i] = row.length
row = row[0]
i += 1
end
# A row vector should be stored as 1xN, not N
#shape.unshift(1) if shape.size == 1
# Then flatten the array.
NMatrix.new(shape, params.flatten, options)
end
#
# call-seq:
# zeros(shape) -> NMatrix
# zeros(shape, dtype: dtype) -> NMatrix
# zeros(shape, dtype: dtype, stype: stype) -> NMatrix
#
# Creates a new matrix of zeros with the dimensions supplied as
# parameters.
#
# * *Arguments* :
# - +shape+ -> Array (or integer for square matrix) specifying the dimensions.
# - +dtype+ -> (optional) Default is +:float64+
# - +stype+ -> (optional) Default is +:dense+.
# * *Returns* :
# - NMatrix filled with zeros.
#
# Examples:
#
# NMatrix.zeros(2) # => 0.0 0.0
# 0.0 0.0
#
# NMatrix.zeros([2, 3], dtype: :int32) # => 0 0 0
# 0 0 0
#
# NMatrix.zeros([1, 5], dtype: :int32) # => 0 0 0 0 0
#
def zeros(shape, opts = {})
NMatrix.new(shape, 0, {:dtype => :float64}.merge(opts))
end
alias :zeroes :zeros
#
# call-seq:
# ones(shape) -> NMatrix
# ones(shape, dtype: dtype, stype: stype) -> NMatrix
#
# Creates a matrix filled with ones.
#
# * *Arguments* :
# - +shape+ -> Array (or integer for square matrix) specifying the shape.
# - +opts+ -> (optional) Hash of options from NMatrix#initialize
# * *Returns* :
# - NMatrix filled with ones.
#
# Examples:
#
# NMatrix.ones([1, 3]) # => 1.0 1.0 1.0
#
# NMatrix.ones([2, 3], dtype: :int32) # => 1 1 1
# 1 1 1
#
def ones(shape, opts={})
NMatrix.new(shape, 1, {:dtype => :float64, :default => 1}.merge(opts))
end
# call-seq:
# ones_like(nm) -> NMatrix
#
# Creates a new matrix of ones with the same dtype and shape as the
# provided matrix.
#
# @param [NMatrix] nm the nmatrix whose dtype and shape will be used
# @return [NMatrix] a new nmatrix filled with ones.
#
def ones_like(nm)
NMatrix.ones(nm.shape, dtype: nm.dtype, stype: nm.stype, capacity: nm.capacity, default: 1)
end
# call-seq:
# zeros_like(nm) -> NMatrix
#
# Creates a new matrix of zeros with the same stype, dtype, and shape
# as the provided matrix.
#
# @param [NMatrix] nm the nmatrix whose stype, dtype, and shape will be used
# @return [NMatrix] a new nmatrix filled with zeros.
#
def zeros_like(nm)
NMatrix.zeros(nm.shape, dtype: nm.dtype, stype: nm.stype, capacity: nm.capacity, default: 0)
end
#
# call-seq:
# eye(shape) -> NMatrix
# eye(shape, dtype: dtype) -> NMatrix
# eye(shape, stype: stype, dtype: dtype) -> NMatrix
#
# Creates an identity matrix (square matrix rank 2).
#
# * *Arguments* :
# - +size+ -> Array (or integer for square matrix) specifying the dimensions.
# - +dtype+ -> (optional) Default is +:float64+
# - +stype+ -> (optional) Default is +:dense+.
# * *Returns* :
# - An identity matrix.
#
# Examples:
#
# NMatrix.eye(3) # => 1.0 0.0 0.0
# 0.0 1.0 0.0
# 0.0 0.0 1.0
#
# NMatrix.eye(3, dtype: :int32) # => 1 0 0
# 0 1 0
# 0 0 1
#
# NMatrix.eye(2, dtype: :int32, stype: :yale) # => 1 0
# 0 1
#
def eye(shape, opts={})
# Fill the diagonal with 1's.
m = NMatrix.zeros(shape, {:dtype => :float64}.merge(opts))
(0...m.shape[0]).each do |i|
m[i, i] = 1
end
m
end
alias :identity :eye
#
# call-seq:
# hilbert(shape) -> NMatrix
# hilbert(shape, dtype: dtype) -> NMatrix
# hilbert(shape, stype: stype, dtype: dtype) -> NMatrix
#
# Creates an hilbert matrix (square matrix).
#
# * *Arguments* :
# - +size+ -> integer ( for square matrix) specifying the dimensions.
# - +dtype+ -> (optional) Default is +:float64+
# - +stype+ -> (optional) Default is +:dense+.
# * *Returns* :
# - A hilbert matrix.
#
# Examples:
#
# NMatrix.hilbert(3) # => 1.0 0.5 0.3333333333333333
# 0.5 0.3333333333333333 0.25
# 0.3333333333333333 0.25 0.2
#
def hilbert(shape, opts={})
m = NMatrix.new([shape,shape], {:dtype => :float64}.merge(opts))
0.upto(shape - 1) do |i|
0.upto(i) do |j|
m[i,j] = 1.0 / (j + i + 1)
m[j,i] = m[i,j] if i != j
end
end
m
end
#
# call-seq:
# inv_hilbert(shape) -> NMatrix
# inv_hilbert(shape, dtype: dtype) -> NMatrix
# inv_hilbert(shape, stype: stype, dtype: dtype) -> NMatrix
#
# Creates an inverse hilbert matrix (square matrix rank 2).
#
# * *Arguments* :
# - +size+ -> Array (or integer for square matrix) specifying the dimensions.
# - +dtype+ -> (optional) Default is +:float64+
# - +stype+ -> (optional) Default is +:dense+.
# * *Returns* :
# - A hilbert matrix.
#
# Examples:
# NMatrix.inv_hilbert(3) # => 9.0, -36.0, 30.0
# -36.0, 192.0, -180.0
# 30.0, -180.0, 180.0
#
#
def inv_hilbert(shape, opts={})
opts = {:dtype => :float64}.merge(opts)
m = NMatrix.new([shape,shape],opts)
combination = NMatrix.new([2*shape,2*shape],opts)
#combinations refers to the combination of n things taken k at a time
0.upto(2*shape-1) do |i|
0.upto(i) do |j|
if j != 0 and j != i
combination[i,j] = combination[i-1,j] + combination[i-1,j-1]
else
combination[i,j] = 1
end
end
end
0.upto(shape-1) do |i|
0.upto(i) do |j|
m[i,j] = combination[shape + j,shape - i - 1] * ((i+j)+1) * \
combination[shape + i,shape - j - 1] * (-1) ** ((i+j)) * \
combination[(i+j),i] * combination[(i+j),i]
m[j,i] = m[i,j] if i != j
end
end
m
end
#
# call-seq:
# diagonals(array) -> NMatrix
# diagonals(array, dtype: dtype, stype: stype) -> NMatrix
#
# Creates a matrix filled with specified diagonals.
#
# * *Arguments* :
# - +entries+ -> Array containing input values for diagonal matrix
# - +options+ -> (optional) Hash with options for NMatrix#initialize
# * *Returns* :
# - NMatrix filled with specified diagonal values.
#
# Examples:
#
# NMatrix.diagonal([1.0,2,3,4]) # => 1.0 0.0 0.0 0.0
# 0.0 2.0 0.0 0.0
# 0.0 0.0 3.0 0.0
# 0.0 0.0 0.0 4.0
#
# NMatrix.diagonal([1,2,3,4], dtype: :int32) # => 1 0 0 0
# 0 2 0 0
# 0 0 3 0
# 0 0 0 4
#
#
def diagonal(entries, opts={})
m = NMatrix.zeros(entries.size,
{:dtype => guess_dtype(entries[0]), :capacity => entries.size + 1}.merge(opts)
)
entries.each_with_index do |n, i|
m[i,i] = n
end
m
end
alias :diag :diagonal
alias :diagonals :diagonal
# Generate a block-diagonal NMatrix from the supplied 2D square matrices.
#
# * *Arguments*
# - +*params+ -> An array that collects all arguments passed to the method. The method
# can receive any number of arguments. Optionally, the last entry of +params+ is
# a hash of options from NMatrix#initialize. All other entries of +params+ are
# the blocks of the desired block-diagonal matrix. Each such matrix block can be
# supplied as a square 2D NMatrix object, or alternatively as an array of arrays
# (with dimensions corresponding to a square matrix), or alternatively as a number.
# * *Returns*
# - NMatrix of block-diagonal form filled with specified matrices
# as the blocks along the diagonal.
#
# * *Example*
#
# a = NMatrix.new([2,2], [1,2,3,4])
# b = NMatrix.new([1,1], [123], dtype: :float64)
# c = Array.new(2) { [[10,10], [10,10]] }
# d = Array[[1,2,3], [4,5,6], [7,8,9]]
# m = NMatrix.block_diagonal(a, b, *c, d, 10.0, 11, dtype: :int64, stype: :yale)
# =>
# [
# [1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
# [3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
# [0, 0, 123, 0, 0, 0, 0, 0, 0, 0, 0, 0]
# [0, 0, 0, 10, 10, 0, 0, 0, 0, 0, 0, 0]
# [0, 0, 0, 10, 10, 0, 0, 0, 0, 0, 0, 0]
# [0, 0, 0, 0, 0, 10, 10, 0, 0, 0, 0, 0]
# [0, 0, 0, 0, 0, 10, 10, 0, 0, 0, 0, 0]
# [0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 0, 0]
# [0, 0, 0, 0, 0, 0, 0, 4, 5, 6, 0, 0]
# [0, 0, 0, 0, 0, 0, 0, 7, 8, 9, 0, 0]
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0]
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11]
# ]
#
def block_diagonal(*params)
options = params.last.is_a?(Hash) ? params.pop : {}
params.each_index do |i|
params[i] = params[i].to_nm if params[i].is_a?(Array) # Convert Array to NMatrix
params[i] = NMatrix.new([1,1], [params[i]]) if params[i].is_a?(Numeric) # Convert number to NMatrix
end
block_sizes = [] #holds the size of each matrix block
params.each do |b|
unless b.is_a?(NMatrix)
raise(ArgumentError, "Only NMatrix or appropriate Array objects or single numbers allowed")
end
raise(ArgumentError, "Only 2D matrices or 2D arrays allowed") unless b.shape.size == 2
raise(ArgumentError, "Only square-shaped blocks allowed") unless b.shape[0] == b.shape[1]
block_sizes << b.shape[0]
end
block_diag_mat = NMatrix.zeros(block_sizes.sum, options)
(0...params.length).each do |n|
# First determine the size and position of the n'th block in the block-diagonal matrix
block_size = block_sizes[n]
block_pos = block_sizes[0...n].sum
# populate the n'th block in the block-diagonal matrix
(0...block_size).each do |i|
(0...block_size).each do |j|
block_diag_mat[block_pos+i,block_pos+j] = params[n][i,j]
end
end
end
return block_diag_mat
end
alias :block_diag :block_diagonal
#
# call-seq:
# random(shape) -> NMatrix
#
# Creates a +:dense+ NMatrix with random numbers between 0 and 1 generated
# by +Random::rand+. The parameter is the dimension of the matrix.
#
# If you use an integer dtype, make sure to specify :scale as a parameter, or you'll
# only get a matrix of 0s.
#
# * *Arguments* :
# - +shape+ -> Array (or integer for square matrix) specifying the dimensions.
# * *Returns* :
# - NMatrix filled with random values.
#
# Examples:
#
# NMatrix.random([2, 2]) # => 0.4859439730644226 0.1783195585012436
# 0.23193766176700592 0.4503345191478729
#
# NMatrix.random([2, 2], :dtype => :byte, :scale => 255) # => [ [252, 108] [44, 12] ]
#
def random(shape, opts={})
scale = opts.delete(:scale) || 1.0
rng = Random.new
random_values = []
# Construct the values of the final matrix based on the dimension.
if opts[:dtype] == :complex64 || opts[:dtype] == :complex128
NMatrix.size(shape).times { |i| random_values << Complex(rng.rand(scale), rng.rand(scale)) }
else
NMatrix.size(shape).times { |i| random_values << rng.rand(scale) }
end
NMatrix.new(shape, random_values, {:dtype => :float64, :stype => :dense}.merge(opts))
end
alias :rand :random
#
# call-seq:
# magic(shape) -> NMatrix
# magic(shape, dtype: dtype) -> NMatrix
#
# The parameter is the dimension of the matrix.
#
# Creates a +:dense+ NMatrix with the following properties:
# - An arrangement of the numbers from 1 to n^2 (n-squared) in the matrix, with each number occurring exactly once.
# - The sum of the entries of any row, any column, or any main diagonal is the same.
# - This sum must be n(n^2+1)/2.
#
# See: http://www.mathworks.com/help/matlab/ref/magic.html
#
# * *Arguments* :
# - +shape+ -> Array (or integer for square matrix) specifying the dimensions.
# - +dtype+ -> (optional) Default is +:float64+
# * *Returns* :
# - NMatrix with the above given properties.
#
# Examples:
#
# NMatrix.magic(3) # => [ [4.0, 9.0, 2.0] [3.0, 5.0, 7.0] [8.0, 1.0, 6.0] ]
#
# NMatrix.magic(4, dtype :int32) # => [ [ 1, 15, 14, 4]
# [12, 6, 7, 9]
# [ 8, 10, 11, 5]
# [13, 3, 2, 16] ]
#
# NMatrix.magic(6,dtype: :int64) # => [ [31, 9, 2, 22, 27, 20]
# [ 3, 32, 7, 21, 23, 25]
# [35, 1, 6, 26, 19, 24]
# [ 4, 36, 29, 13, 18, 11]
# [30, 5, 34, 12, 14, 16]
# [ 8, 28, 33, 17, 10, 15] ]
#
def magic(shape, opts={})
raise(ArgumentError, "shape of two is not allowed") if shape == 2
nm = NMatrix.new([shape,shape], 0, {:dtype => :float64}.merge(opts))
if shape % 2 != 0
MagicHelpers.odd_magic nm, shape
elsif shape % 4 == 0
MagicHelpers.doubly_even_magic nm, shape
else
MagicHelpers.singly_even_magic nm, shape
end
nm
end
#
# call-seq:
# linspace(base, limit) -> 1x100 NMatrix
# linspace(base, limit, *shape) -> NMatrix
#
# Returns an NMatrix with +[shape[0] x shape[1] x .. x shape[dim-1]]+ values of dtype +:float64+ equally spaced from
# +base+ to +limit+, inclusive.
#
# See: http://www.mathworks.com/help/matlab/ref/linspace.html
#
# * *Arguments* :
# - +base+ -> The first value in the sequence.
# - +limit+ -> The last value in the sequence.
# - +shape+ -> Desired output shape. Default returns a 1x100 row vector.
# * *Returns* :
# - NMatrix with +:float64+ values.
#
# Examples :-
#
# NMatrix.linspace(1,Math::PI, 6)
# =>[1.0,
# 1.4283185005187988,
# 1.8566370010375977,
# 2.2849555015563965,
# 2.7132740020751953,
# 3.1415927410125732
# ]
#
# NMatrix.linspace(1,10, [3,2])
# =>[
# [ 1.0, 2.799999952316284]
# [4.599999904632568, 6.400000095367432]
# [8.199999809265137, 10.0]
# ]
#
def linspace(base, limit, shape = [100])
# Convert shape to array format
shape = [shape] if shape.is_a? Integer
#Calculate number of elements
count = shape.inject(:*)
# Linear spacing between elements calculated in step
# step = limit - base / (count - 1)
# [Result Sequence] = [0->N sequence] * step + [Base]
step = (limit - base) * (1.0 / (count - 1))
result = NMatrix.seq(shape, {:dtype => :float64}) * step
result += NMatrix.new(shape, base)
result
end
# call-seq:
# logspace(base, limit) -> 1x50 NMatrix with exponent_base = 10
# logspace(base, limit, shape , exponent_base:) -> NMatrix
# logspace(base, :pi, n) -> 1xn NMatrix with interval [10 ^ base, Math::PI]
#
# Returns an NMatrix with +[shape[0] x shape[1] x .. x shape[dim-1]]+ values of dtype +:float64+ logarithmically spaced from
# +exponent_base ^ base+ to +exponent_base ^ limit+, inclusive.
#
# See: http://www.mathworks.com/help/matlab/ref/logspace.html
#
# * *Arguments* :
# - +base+ -> exponent_base ** base is the first value in the sequence
# - +limit+ -> exponent_base ** limit is the last value in the sequence.
# - +shape+ -> Desired output shape. Default returns a 1x50 row vector.
# * *Returns* :
# - NMatrix with +:float64+ values.
#
# Examples :-
#
# NMatrix.logspace(1,:pi,7)
# =>[
# 10.0000,
# 8.2450,
# 6.7980,
# 5.6050,
# 4.6213,
# 3.8103,
# 3.1416
# ]
#
# NMatrix.logspace(1,2,[3,2])
# =>[
# [10.0, 15.8489]
# [25.1189, 39.8107]
# [63.0957, 100.0]
# ]
#
def logspace(base, limit, shape = [50], exponent_base: 10)
#Calculate limit for [10 ^ base ... Math::PI] if limit = :pi
limit = Math.log(Math::PI, exponent_base = 10) if limit == :pi
shape = [shape] if shape.is_a? Integer
#[base...limit] -> [exponent_base ** base ... exponent_base ** limit]
result = NMatrix.linspace(base, limit, shape)
result.map {|element| exponent_base ** element}
end
#
# call-seq:
# linspace(base, limit) -> 1x100 NMatrix
# linspace(base, limit, *shape) -> NMatrix
#
# Returns an NMatrix with +[shape[0] x shape[1] x .. x shape[dim-1]]+ values of dtype +:float64+ equally spaced from
# +base+ to +limit+, inclusive.
#
# See: http://www.mathworks.com/help/matlab/ref/linspace.html
#
# * *Arguments* :
# - +base+ -> The first value in the sequence.
# - +limit+ -> The last value in the sequence.
# - +shape+ -> Desired output shape. Default returns a 1x100 row vector.
# * *Returns* :
# - NMatrix with +:float64+ values.
#
# Examples :-
#
# NMatrix.linspace(1,Math::PI, 6)
# =>[1.0,
# 1.4283185005187988,
# 1.8566370010375977,
# 2.2849555015563965,
# 2.7132740020751953,
# 3.1415927410125732
# ]
#
# NMatrix.linspace(1,10, [3,2])
# =>[
# [ 1.0, 2.799999952316284]
# [4.599999904632568, 6.400000095367432]
# [8.199999809265137, 10.0]
# ]
#
def linspace(base, limit, shape = [100])
# Convert shape to array format
shape = [shape] if shape.is_a? Integer
#Calculate number of elements
count = shape.inject(:*)
# Linear spacing between elements calculated in step
# step = limit - base / (count - 1)
# [Result Sequence] = [0->N sequence] * step + [Base]
step = (limit - base) * (1.0 / (count - 1))
result = NMatrix.seq(shape, {:dtype => :float64}) * step
result += NMatrix.new(shape, base)
result
end
# call-seq:
# logspace(base, limit) -> 1x50 NMatrix with exponent_base = 10
# logspace(base, limit, shape , exponent_base:) -> NMatrix
# logspace(base, :pi, n) -> 1xn NMatrix with interval [10 ^ base, Math::PI]
#
# Returns an NMatrix with +[shape[0] x shape[1] x .. x shape[dim-1]]+ values of dtype +:float64+ logarithmically spaced from
# +exponent_base ^ base+ to +exponent_base ^ limit+, inclusive.
#
# See: http://www.mathworks.com/help/matlab/ref/logspace.html
#
# * *Arguments* :
# - +base+ -> exponent_base ** base is the first value in the sequence
# - +limit+ -> exponent_base ** limit is the last value in the sequence.
# - +shape+ -> Desired output shape. Default returns a 1x50 row vector.
# * *Returns* :
# - NMatrix with +:float64+ values.
#
# Examples :-
#
# NMatrix.logspace(1,:pi,7)
# =>[
# 10.0000,
# 8.2450,
# 6.7980,
# 5.6050,
# 4.6213,
# 3.8103,
# 3.1416
# ]
#
# NMatrix.logspace(1,2,[3,2])
# =>[
# [10.0, 15.8489]
# [25.1189, 39.8107]
# [63.0957, 100.0]
# ]
#
def logspace(base, limit, shape = [50], exponent_base: 10)
#Calculate limit for [10 ^ base ... Math::PI] if limit = :pi
limit = Math.log(Math::PI, exponent_base = 10) if limit == :pi
shape = [shape] if shape.is_a? Integer
#[base...limit] -> [exponent_base ** base ... exponent_base ** limit]
result = NMatrix.linspace(base, limit, shape)
result.map {|element| exponent_base ** element}
end
#
# call-seq:
# seq(shape) -> NMatrix
# seq(shape, options) -> NMatrix
# bindgen(shape) -> NMatrix of :byte
# indgen(shape) -> NMatrix of :int64
# findgen(shape) -> NMatrix of :float32
# dindgen(shape) -> NMatrix of :float64
# cindgen(shape) -> NMatrix of :complex64
# zindgen(shape) -> NMatrix of :complex128
# rbindgen(shape) -> NMatrix of :object
#
# Creates a matrix filled with a sequence of integers starting at zero.
#
# * *Arguments* :
# - +shape+ -> Array (or integer for square matrix) specifying the dimensions.
# - +options+ -> (optional) Options permissible for NMatrix#initialize
# * *Returns* :
# - NMatrix filled with values 0 through +size+.
#
# Examples:
#
# NMatrix.seq(2) # => 0 1
# 2 3
#
# NMatrix.seq([3, 3], dtype: :float32) # => 0.0 1.0 2.0
# 3.0 4.0 5.0
# 6.0 7.0 8.0
#
def seq(shape, options={})
# Construct the values of the final matrix based on the dimension.
values = (0 ... NMatrix.size(shape)).to_a
# It'll produce :int32, except if a dtype is provided.
NMatrix.new(shape, values, {:stype => :dense}.merge(options))
end
{:bindgen => :byte, :indgen => :int64, :findgen => :float32, :dindgen => :float64,
:cindgen => :complex64, :zindgen => :complex128,
:rbindgen => :object}.each_pair do |meth, dtype|
define_method(meth) { |shape| NMatrix.seq(shape, :dtype => dtype) }
end
end
end
module NVector #:nodoc:
class << self
#
# call-seq:
# new(shape) -> NVector
# new(stype, shape) -> NVector
# new(shape, init) -> NVector
# new(:dense, shape, init) -> NVector
# new(:list, shape, init) -> NVector
# new(shape, init, dtype) -> NVector
# new(stype, shape, init, dtype) -> NVector
# new(stype, shape, dtype) -> NVector
#
# Creates a new NVector. See also NMatrix#initialize for a more detailed explanation of
# the arguments.
#
# * *Arguments* :
# - +stype+ -> (optional) Storage type of the vector (:list, :dense, :yale). Defaults to :dense.
# - +shape+ -> Shape of the vector. Accepts [n,1], [1,n], or n, where n is a Fixnum.
# - +init+ -> (optional) Yale: capacity; List: default value (0); Dense: initial value or values (uninitialized by default).
# - +dtype+ -> (optional if +init+ provided) Data type stored in the vector. For :dense and :list, can be inferred from +init+.
# * *Returns* :
# -
#
def new(*args)
stype = args[0].is_a?(Symbol) ? args.shift : :dense
shape = args[0].is_a?(Array) ? args.shift : [1,args.shift]
if shape.size != 2 || !shape.include?(1) || shape == [1,1]
raise(ArgumentError, "shape must be a Fixnum or an Array of positive Fixnums where exactly one value is 1")
end
warn "NVector is deprecated and not guaranteed to work any longer"
NMatrix.new(stype, shape, *args)
end
#
# call-seq:
# zeros(size) -> NMatrix
# zeros(size, dtype) -> NMatrix
#
# Creates a new vector of zeros with the dimensions supplied as
# parameters.
#
# * *Arguments* :
# - +size+ -> Array (or integer for square matrix) specifying the dimensions.
# - +dtype+ -> (optional) Default is +:float64+.
# * *Returns* :
# - NVector filled with zeros.
#
# Examples:
#
# NVector.zeros(2) # => 0.0
# 0.0
#
# NVector.zeros(3, :int32) # => 0
# 0
# 0
#
def zeros(size, dtype = :float64)
NMatrix.new([size,1], 0, dtype: dtype)
end
alias :zeroes :zeros
#
# call-seq:
# ones(size) -> NVector
# ones(size, dtype) -> NVector
#
# Creates a vector of ones with the dimensions supplied as
# parameters.
#
# * *Arguments* :
# - +size+ -> Array (or integer for square matrix) specifying the dimensions.
# - +dtype+ -> (optional) Default is +:float64+.
# * *Returns* :
# - NVector filled with ones.
#
# Examples:
#
# NVector.ones(2) # => 1.0
# 1.0
#
# NVector.ones(3, :int32) # => 1
# 1
# 1
#
def ones(size, dtype = :float64)
NMatrix.new([size,1], 1, dtype: dtype)
end
#
# call-seq:
# random(size) -> NVector
#
# Creates a vector with random numbers between 0 and 1 generated by
# +Random::rand+ with the dimensions supplied as parameters.
#
# * *Arguments* :
# - +size+ -> Array (or integer for square matrix) specifying the dimensions.
# - +opts+ -> (optional) NMatrix#initialize options
# * *Returns* :
# - NVector filled with random numbers generated by the +Random+ class.
#
# Examples:
#
# NVector.rand(2) # => 0.4859439730644226
# 0.1783195585012436
#
def random(size, opts = {})
rng = Random.new
random_values = []
size.times { |i| random_values << rng.rand }
NMatrix.new([size,1], random_values, opts)
end
#
# call-seq:
# seq(n) -> NVector
# seq(n, dtype) -> NVector
#
# Creates a vector with a sequence of +n+ integers starting at zero. You
# can choose other types based on the dtype parameter.
#
# * *Arguments* :
# - +n+ -> Number of integers in the sequence.
# - +dtype+ -> (optional) Default is +:int64+.
# * *Returns* :
# - NVector filled with +n+ integers.
#
# Examples:
#
# NVector.seq(2) # => 0
# 1
#
# NVector.seq(3, :float32) # => 0.0
# 1.0
# 2.0
#
def seq(size, dtype = :int64)
values = (0 ... size).to_a
NMatrix.new([size,1], values, dtype: dtype)
end
#
# call-seq:
# indgen(n) -> NVector
#
# Returns an integer NVector. Equivalent to seq(n, :int32).
#
# * *Arguments* :
# - +n+ -> Size of the sequence.
# * *Returns* :
# - NVector filled with +n+ integers of dtype +:int32+.
#
def indgen(n)
NVector.seq(n, :int32)
end
#
# call-seq:
# findgen(n) -> NVector
#
# Returns a float NVector. Equivalent to seq(n, :float32).
#
# * *Arguments* :
# - +n+ -> Size of the sequence.
# * *Returns* :
# - NVector filled with +n+ integers of dtype +:float32+.
#
def findgen(n)
NVector.seq(n, :float32)
end
#
# call-seq:
# bindgen(n) -> NVector
#
# Returns a byte NVector. Equivalent to seq(n, :byte).
#
# * *Arguments* :
# - +n+ -> Size of the sequence.
# * *Returns* :
# - NVector filled with +n+ integers of dtype +:byte+.
#
def bindgen(n)
NVector.seq(n, :byte)
end
#
# call-seq:
# cindgen(n) -> NVector
#
# Returns a complex NVector. Equivalent to seq(n, :complex64).
#
# * *Arguments* :
# - +n+ -> Size of the sequence.
# * *Returns* :
# - NVector filled with +n+ integers of dtype +:complex64+.
#
def cindgen(n)
NVector.seq(n, :complex64)
end
#
# call-seq:
# linspace(a, b) -> NVector
# linspace(a, b, n) -> NVector
#
# Returns a NVector with +n+ values of dtype +:float64+ equally spaced from
# +a+ to +b+, inclusive.
#
# See: http://www.mathworks.com/help/matlab/ref/linspace.html
#
# * *Arguments* :
# - +a+ -> The first value in the sequence.
# - +b+ -> The last value in the sequence.
# - +n+ -> The number of elements. Default is 100.
# * *Returns* :
# - NVector with +n+ +:float64+ values.
#
# Example:
# x = NVector.linspace(0, Math::PI, 1000)
# x.pretty_print
# [0.0
# 0.0031447373909807737
# 0.006289474781961547
# ...
# 3.135303178807831
# 3.138447916198812
# 3.141592653589793]
# => nil
#
def linspace(a, b, n = 100)
# Formula: seq(n) * step + a
# step = ((b - a) / (n - 1))
step = (b - a) * (1.0 / (n - 1))
# dtype = :float64 is used to prevent integer coercion.
result = NVector.seq(n, :float64) * NMatrix.new([n,1], step, dtype: :float64)
result += NMatrix.new([n,1], a, dtype: :float64)
result
end
#
# call-seq:
# logspace(a, b) -> NVector
# logspace(a, b, n) -> NVector
#
# Returns a NVector with +n+ values of dtype +:float64+ logarithmically
# spaced from +10^a+ to +10^b+, inclusive.
#
# See: http://www.mathworks.com/help/matlab/ref/logspace.html
#
# * *Arguments* :
# - +a+ -> The first value in the sequence.
# - +b+ -> The last value in the sequence.
# - +n+ -> The number of elements. Default is 100.
# * *Returns* :
# - NVector with +n+ +:float64+ values.
#
# Example:
# x = NVector.logspace(0, Math::PI, 10)
# x.pretty_print
# [1.0
# 2.2339109164570266
# 4.990357982665873
# 11.148015174505757
# 24.903672795156997
# 55.632586516975095
# 124.27824233101062
# 277.6265222213364
# 620.1929186882427
# 1385.4557313670107]
# => nil
#
def logspace(a, b, n = 100)
# Formula: 10^a, 10^(a + step), ..., 10^b, where step = ((b-a) / (n-1)).
result = NVector.linspace(a, b, n)
result.each_stored_with_index { |element, i| result[i] = 10 ** element }
result
end
end
end
# This constant is intended as a simple constructor for NMatrix meant for
# experimenting.
#
# Examples:
#
# a = N[ 1,2,3,4 ] => 1 2 3 4
#
# a = N[ 1,2,3,4, :int32 ] => 1 2 3 4
#
# a = N[ [1,2,3], [3,4,5] ] => 1 2 3
# 3 4 5
#
# a = N[ 3,6,9 ].transpose => 3
# 6
# 9
N = NMatrix