Class: NMatrix

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

Overview

require ‘pp’

Instance Method Summary collapse

Instance Method Details

#cosine_normalization ⇒ Object 



172
173
174
175
176
177
178
179
180
181
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 172

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



9
10
11
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 9

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 



13
14
15
16
17
18
19
20
21
22
23
24
25
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 13

def div_by_vector(vector, by=:col)
	new_matrix =  NMatrix.zeros(self.shape, dtype: self.dtype)
	if by == :col
		self.cols.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 



80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 80

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 



27
28
29
30
31
32
33
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 27

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-12) ⇒ Object

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



56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 56

def max_eigenvalue(n=100, error = 10e-12) # do not set error too low or the eigenvalue cannot stabilised around the real one
	max_eigenvalue = 0.0
	length = self.cols
	v = NMatrix.random([self.cols, 1], dtype: self.dtype)
	# http://web.mit.edu/18.06/www/Spring17/Power-Method.pdf 
	#IMPLEMENTATION PROBLEM: RESULTS ARE TOO VARIABLE
	last_max_eigenvalue = nil
	n.times do
		v = self.dot(v) # calculate the matrix-by-vector product Mv
		v = v / v.frobenius_norm # calculate the norm and normalize the vector
		max_eigenvalue = v.vector_product(self.dot(v)) / v.vector_self_product #Rayleigh quotient 
		# Rayleigh quotient: lambda = vMv/vv
		# v is a vector so vv is inner product of one vector with self (use vector_self_product); 
		# Mv gives a vector, so vMv is the inner product of two different vectors (use vector_product)
		break if !last_max_eigenvalue.nil? && last_max_eigenvalue - max_eigenvalue <= 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



35
36
37
38
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 35

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 



76
77
78
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 76

def min_eigenvalue(n=100, error = 10e-12)
	return  self.invert.max_eigenvalue(n, error)
end

#vector_product(vec_b) ⇒ Object 



40
41
42
43
44
45
46
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 40

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 



48
49
50
51
52
53
54
 
# File 'lib/expcalc/nmatrix_expansion.rb', line 48

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