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



102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
 
# File 'lib/expcalc/numo_expansion.rb', line 102

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

#bmatrix_rectangular_to_hash(x_names, y_names) ⇒ Object 



147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
 
# File 'lib/expcalc/numo_expansion.rb', line 147

def bmatrix_rectangular_to_hash(x_names, y_names)
	hash = {}
	x_names.each_with_index do |x_name, x_pos|
		y_names.each_with_index do |y_name, y_pos|
			associationValue = self[x_pos, y_pos]
			if associationValue > 0
				query = hash[x_name]
				if query.nil?
					hash[x_name] = [y_name]
				else
					hash[x_name] << y_name
				end
				query = hash[y_name]
				if query.nil?
					hash[y_name] = [x_name]
				else
					hash[y_name] << x_name
				end
			end
		end
	end
	return hash
end

#bmatrix_squared_to_hash(elements_names) ⇒ Object 



129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
 
# File 'lib/expcalc/numo_expansion.rb', line 129

def bmatrix_squared_to_hash(elements_names)
	hash = {}
	elements_names.each_with_index do |x_name, x_pos|
		elements_names.each_with_index do |y_name, y_pos|
			associationValue = self[x_pos, y_pos]
			if associationValue > 0
				query = hash[x_name]
				if query.nil?
					hash[x_name] = [y_name]
				else
					hash[x_name] << y_name
				end
			end
		end
	end
	return hash
end

#cosine_normalization ⇒ Object 



372
373
374
375
376
377
378
379
380
 
# File 'lib/expcalc/numo_expansion.rb', line 372

def cosine_normalization
	normalized_matrix =  Numo::DFloat.zeros(self.shape) 
	self.each_with_index 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



213
214
215
 
# File 'lib/expcalc/numo_expansion.rb', line 213

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 



217
218
219
220
221
222
223
224
225
226
227
228
229
 
# File 'lib/expcalc/numo_expansion.rb', line 217

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 



280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
 
# File 'lib/expcalc/numo_expansion.rb', line 280

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 



231
232
233
234
235
236
237
 
# File 'lib/expcalc/numo_expansion.rb', line 231

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



260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
 
# File 'lib/expcalc/numo_expansion.rb', line 260

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



239
240
241
242
 
# File 'lib/expcalc/numo_expansion.rb', line 239

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 



276
277
278
 
# File 'lib/expcalc/numo_expansion.rb', line 276

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 



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

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 



244
245
246
247
248
249
250
 
# File 'lib/expcalc/numo_expansion.rb', line 244

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 



252
253
254
255
256
257
258
 
# File 'lib/expcalc/numo_expansion.rb', line 252

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

#wmatrix_rectangular_to_hash(x_names, y_names) ⇒ Object 



189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
 
# File 'lib/expcalc/numo_expansion.rb', line 189

def wmatrix_rectangular_to_hash(x_names, y_names)
	hash = {}
	x_names.each_with_index do |x_name, x_pos|
		y_names.each_with_index do |y_name, y_pos|
			associationValue = self[x_pos, y_pos].to_f
			if associationValue > 0
				query = hash[x_name]
				if query.nil?
					hash[x_name] = {y_name => associationValue}
				else
					hash[x_name][y_name] = associationValue
				end
				query = hash[y_name]
				if query.nil?
					hash[y_name] = {x_name => associationValue}
				else
					hash[y_name][x_name] = associationValue
				end
			end
		end
	end
	return hash
end

#wmatrix_squared_to_hash(layers) ⇒ Object 



171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
 
# File 'lib/expcalc/numo_expansion.rb', line 171

def wmatrix_squared_to_hash(layers)
	hash = {}
	layers.each_with_index do |x_name, x_pos|
		layers.each_with_index do |y_name, y_pos|
			associationValue = self[x_pos, y_pos]
			if associationValue > 0
				query = hash[x_name]
				if query.nil?
					hash[x_name] = {y_name => associationValue}
				else
					hash[x_name][y_name] = associationValue
				end
			end
		end
	end
	return hash
end