Class: Numo::NArray

Inherits:
Object
  • Object
show all
Defined in:
lib/expcalc/numo_expansion.rb

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.load(input_file, type = 'txt', splitChar = "\t") ⇒ Object

TODO: load and save must have exact coherence between load ouput and save input Take into acount github.com/ankane/npy to load and read mats it allows serialize matrices larger than ruby marshal format



71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
 
# File 'lib/expcalc/numo_expansion.rb', line 71

def self.load(input_file, type='txt', splitChar="\t")
	matrix = nil
	if type == 'txt' 
		counter = 0
		File.open(input_file).each do |line|
	    	line.chomp!
    		row = line.split(splitChar).map{|c| c.to_f }
    		if matrix.nil?
    			matrix = Numo::DFloat.zeros(row.length, row.length)
    		end
    		row.each_with_index do |val, i|
    			matrix[counter, i] = val if val != 0
    		end
    		counter += 1
		end
	elsif type == 'npy'
		matrix = Npy.load(input_file)
	end
	return matrix
end

Instance Method Details

#cosine_normalization ⇒ Object 



258
259
260
261
262
263
264
265
266
267
 
# File 'lib/expcalc/numo_expansion.rb', line 258

def cosine_normalization
	normalized_matrix =  NMatrix.zeros(self.shape, dtype: self.dtype)
	#normalized_matrix =  NMatrix.zeros(self.shape, dtype: :complex64)
	self.each_with_indices do |val, i, j|
		norm = val/CMath.sqrt(self[i, i] * self[j,j])
		#abort("#{norm} has non zero imaginary part" ) if norm.imag != 0
		normalized_matrix[i, j] = norm#.real
	end
	return normalized_matrix
end

#div(second_mat) ⇒ Object

Matrix division A/B => A.dot(B.pinv) #stackoverflow.com/questions/49225693/matlab-matrix-division-into-python



99
100
101
 
# File 'lib/expcalc/numo_expansion.rb', line 99

def div(second_mat) #Matrix division A/B => A.dot(B.pinv) #https://stackoverflow.com/questions/49225693/matlab-matrix-division-into-python
	return self.dot(second_mat.pinv)
end

#div_by_vector(vector, by = :col) ⇒ Object 



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

def div_by_vector(vector, by=:col)
	new_matrix =  self.new_zeros
	if by == :col
		self.shape.last.times do |n|
			vector.each_with_indices do |val, i, j|
				new_matrix[i, n] = self[i, n].fdiv(val)
			end
		end
	elsif by == :row

	end
	return new_matrix
end

#expm ⇒ Object 



166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
 
# File 'lib/expcalc/numo_expansion.rb', line 166

def expm
	return compute_py_method{|mat| expm(mat)}
	#return compute_py_method(self){|mat| expm(mat)}
	##################################################
	# matlab pade aproximation 
	################################################
	### toolbox/matlab/demos/expmdemo1.m (Golub and Van Loan, Matrix Computations, Algorithm 11.3-1.)		
	
	#fraction, exponent = Math.frexp(max_norm)
	#s = [0, exponent+1].max
	#a = self/2**s

	## Pade approximation for exp(A)
	#x = a
	#c = 0.5
	#ac = a*c
	#e = NMatrix.identity(a.shape, dtype: a.dtype) + ac
	#d = NMatrix.identity(a.shape, dtype: a.dtype) - ac
	#q = 6
	#p = true
	#(2..q).each do |k|
	#	c = c * (q-k+1) / (k*(2*q-k+1))
	#	x = a.dot(x)
	#	cX =  x * c
	#	e = e + cX
	#	if p
	#		d = d + cX
	#	else
	#		d = d - cX
	#	end
	#	p = !p
	#end
	#e = d.solve(e) #solve

	## Undo scaling by repeated squaring
	#(1..s).each do
	#	e = e.dot(e) 
	#end
	#return e

	###################################
	## Old python Pade aproximation
	###################################
	#### Pade aproximation: https://github.com/rngantner/Pade_PyCpp/blob/master/src/expm.py
	#a_l1 = max_norm
	#n_squarings = 0
	#if self.dtype == :float64 || self.dtype == :complex128
	#	if a_l1 < 1.495585217958292e-002
	#		u,v = _pade3(self)
        #elsif a_l1 < 2.539398330063230e-001
	#		u,v = _pade5(self)
        #elsif a_l1 < 9.504178996162932e-001
	#		u,v = _pade7(self)
        #elsif a_l1 < 2.097847961257068e+000
	#		u,v = _pade9(self)
	#	else
	#		maxnorm = 5.371920351148152
	#		n_squarings = [0, Math.log2(a_l1 / maxnorm).ceil].max
	#		mat = self / 2**n_squarings
	#		u,v = _pade13(mat)
	#	end
	#elsif self.dtype == :float32 || self.dtype == :complex64
	#	if a_l1 < 4.258730016922831e-001
	#		u,v = _pade3(self)
	#    elsif a_l1 < 1.880152677804762e+000
	#		u,v = _pade5(self)
	#	else
	#		maxnorm = 3.925724783138660
	#		n_squarings = [0, Math.log2(a_l1 / maxnorm).ceil].max
	#		mat = self / 2**n_squarings
	#		u,v = _pade7(mat)
	#	end
	#end
	#p = u + v
	#q = -u + v
	#r = q.solve(p)
	#n_squarings.times do
	#	r = r.dot(r)
	#end
	#return r

	######################
	# Exact computing
	######################
	#####expm(matrix) = V*diag(exp(diag(D)))/V; V => eigenvectors(right), D => eigenvalues (right). # https://es.mathworks.com/help/matlab/ref/expm.html
	#eigenvalues, eigenvectors = NMatrix::LAPACK.geev(self, :right)
	#eigenvalues.map!{|val| Math.exp(val)}
	#numerator = eigenvectors.dot(NMatrix.diagonal(eigenvalues, dtype: self.dtype))
	#matrix_exp = numerator.div(eigenvectors)
	#return matrix_exp
end

#frobenius_norm ⇒ Object 



117
118
119
120
121
122
123
 
# File 'lib/expcalc/numo_expansion.rb', line 117

def frobenius_norm
	fro = 0.0
	self.each do |value|
		fro += value.abs ** 2
	end
	return fro ** 0.5
end

#max_eigenvalue(n = 100, error = 10e-15) ⇒ Object

do not set error too low or the eigenvalue cannot stabilised around the real one



146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
 
# File 'lib/expcalc/numo_expansion.rb', line 146

def max_eigenvalue(n=100, error = 10e-15) # do not set error too low or the eigenvalue cannot stabilised around the real one
	max_eigenvalue = 0.0
	length = self.shape.last
	v = Numo::DFloat.new(length).rand
	# http://web.mit.edu/18.06/www/Spring17/Power-Method.pdf 
	last_max_eigenvalue = nil
	n.times do
		v = Numo::Linalg.dot(self, v) 
		v = v / Numo::Linalg.norm(v) 
		max_eigenvalue = Numo::Linalg.dot(v, Numo::Linalg.dot(self, v)) / Numo::Linalg.dot(v,v) #Rayleigh quotient
		break if !last_max_eigenvalue.nil? && (last_max_eigenvalue - max_eigenvalue).abs <= error
		last_max_eigenvalue = max_eigenvalue
	end
	return max_eigenvalue
end

#max_norm ⇒ Object

docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.norm.html, ord parameter = 1



125
126
127
128
 
# File 'lib/expcalc/numo_expansion.rb', line 125

def  max_norm #https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.norm.html, ord parameter = 1
	sums = self.abs.sum(1)
	return sums.max[0, 0]
end

#min_eigenvalue(n = 100, error = 10e-12) ⇒ Object 



162
163
164
 
# File 'lib/expcalc/numo_expansion.rb', line 162

def min_eigenvalue(n=100, error = 10e-12)
	return  Numo::Linalg.inv(self).max_eigenvalue(n, error)
end

#save(matrix_filename, x_axis_names = nil, x_axis_file = nil, y_axis_names = nil, y_axis_file = nil) ⇒ Object 



92
93
94
95
96
 
# File 'lib/expcalc/numo_expansion.rb', line 92

def save(matrix_filename, x_axis_names=nil, x_axis_file=nil, y_axis_names=nil, y_axis_file=nil)
  File.open(x_axis_file, 'w'){|f| f.print x_axis_names.join("\n") } if !x_axis_names.nil?
  File.open(y_axis_file, 'w'){|f| f.print y_axis_names.join("\n") } if !y_axis_names.nil?
  Npy.save(matrix_filename, self)
end

#vector_product(vec_b) ⇒ Object 



130
131
132
133
134
135
136
 
# File 'lib/expcalc/numo_expansion.rb', line 130

def vector_product(vec_b)
	product = 0.0
	self.each_with_indices do |val, i, j|
		product += val * vec_b[i, j]
	end
	return product
end

#vector_self_product ⇒ Object 



138
139
140
141
142
143
144
 
# File 'lib/expcalc/numo_expansion.rb', line 138

def vector_self_product
	product = 0.0
	self.each_stored_with_indices do |val, i, j|
		product +=  val ** 2
	end
	return product
end