ssa

1
package oss
2

3
/*
4
Пакет для работы с матрицами.
5
ПОдстоянно должен дорабатываться, в соответствии с gonum
6
*/
7

8
import (
9
	"errors"
10
	"fmt"
11
	"math"
12
	"sort"
13
	"strings"
14

15
	"gonum.org/v1/gonum/mat"
16
)
17

18
// Горизонтальная сумма массива
19
func Sum2(a mat.Dense) []float64 {
20
	var output []float64
21
	r, _ := a.Dims()
22
	for i := 0; i < r; i++ {
23
		output = append(output, mat.Sum(a.RowView(i)))
24
	}
25
	return output
26
}
27

28
// МАксимальный элемент массива
29
func Max(arr []float64) float64 {
30
	if len(arr) == 0 {
31
		return 0.0
32
	}
33
	max_num := arr[0]
34
	for i := 0; i < len(arr); i++ {
35
		if arr[i] > max_num {
36
			max_num = arr[i]
37
		}
38
	}
39
	return max_num
40
}
41

42
func ATa(matr mat.Dense) mat.Dense { // Multipy matrix AT*A
43
	a := mat.Matrix(&matr)
44
	aT := a.T()
45
	ad := mat.DenseCopyOf(a)
46
	aTd := mat.DenseCopyOf(aT)
47
	n1, _ := aTd.Dims()
48
	_, m2 := ad.Dims()
49
	output := mat.NewDense(n1, m2, nil)
50
	output.Mul(aTd, ad)
51
	return *output
52
}
53

54
func MulVecToVec(a, b mat.Vector) float64 {
55
	a_matrix := mat.Matrix(a).T()
56
	var varVecDense mat.VecDense
57
	varVecDense.MulVec(a_matrix, b)
58
	//fmt.Println(varVecDense)
59
	return varVecDense.AtVec(0)
60
}
61

62
// модуль от всех значений вектора
63
func AbsVector(vect mat.VecDense) mat.VecDense {
64
	for i := 0; i < vect.Len(); i++ {
65
		vect.SetVec(i, math.Abs(vect.AtVec(i)))
66
	}
67
	return vect
68
}
69

70
/*
71

72
// Среднее значение каждого столбца
73
func mean(m mat.Dense) []float64 {
74
	_, c_m := m.Dims()
75
	outputArray := make([]float64, c_m)
76
	for ind := range outputArray {
77
		vect := colDense(m, ind)
78
		outputArray[ind] = averge(vect)
79
	}
80
	return outputArray
81
}
82
*/
83

84
// Среднее массива float64
85
func Averge(array []float64) float64 {
86
	var sum float64
87
	for _, val := range array {
88
		sum += val
89
	}
90
	return sum / float64(len(array))
91
}
92

93
// Вернуть subdiagonal
94
func Subdiagonal(m mat.Dense, k int) ([]float64, error) {
95
	var outputArray []float64
96
	r_m, c_m := m.Dims()
97
	if k == 0 {
98
		for i := 0; i < r_m && i < c_m; i++ {
99
			outputArray = append(outputArray, m.At(i, i))
100
		}
101
	} else if k < 0 {
102
		if k > r_m {
103
			return nil, errors.New("k > matrix ROW")
104
		}
105
		for i := 0; i-k < r_m && i < c_m; i++ {
106
			outputArray = append(outputArray, m.At(i-k, i))
107
		}
108
	} else if k > 0 {
109
		if k > c_m {
110
			return nil, errors.New("k > matrix COL")
111
		}
112
		for i := 0; i+k < c_m && i < r_m; i++ {
113
			outputArray = append(outputArray, m.At(i, i+k))
114
		}
115
	}
116
	return outputArray, nil
117
}
118

119
// Вернуть колонку из матрицы
120
func ColDense(m mat.Dense, ind int) []float64 {
121
	row, _ := m.Dims()
122
	outputArray := make([]float64, row)
123
	for i := 0; i < row; i++ {
124
		outputArray[i] = m.At(i, ind)
125
	}
126
	return outputArray
127
}
128

129
// Вернуть строку из матрицы
130
func RowDense(m mat.Dense, ind int) []float64 {
131
	_, col := m.Dims()
132
	outputArray := make([]float64, col)
133
	for j := 0; j < col; j++ {
134
		outputArray[j] = m.At(ind, j)
135
	}
136
	return outputArray
137
}
138

139
func Diag(mat mat.Dense, R int) []float64 {
140
	ret := make([]float64, R)
141
	for ind := range ret {
142
		ret[ind] = mat.At(ind, ind)
143
	}
144
	return ret
145
}
146

147
func Make_diag_danse(arr []float64) mat.Dense {
148
	lensOfArray := len(arr)
149
	dens := mat.NewDense(lensOfArray, lensOfArray, nil)
150
	for i := 0; i < len(arr); i++ {
151
		dens.Set(i, i, arr[i])
152
	}
153
	return *dens
154
}
155

156
func VecDense_in_float64(vec mat.VecDense) []float64 {
157
	leng, _ := vec.Dims()
158
	output := make([]float64, leng)
159
	for ind := range output {
160
		output[ind] = vec.AtVec(ind)
161
	}
162
	return output
163
}
164
func Vec_in_float64(vec mat.Vector) []float64 {
165
	leng, _ := vec.Dims()
166
	output := make([]float64, leng)
167
	for ind := range output {
168
		output[ind] = vec.AtVec(ind)
169
	}
170
	return output
171
}
172

173
// AsSymDense attempts return a SymDense from the provided Dense.
174
func AsSymDense(m *mat.Dense) (*mat.SymDense, error) {
175
	r, c := m.Dims()
176
	if r != c {
177
		return nil, errors.New("matrix must be square")
178
	}
179
	mT := m.T()
180
	vals := make([]float64, r*c)
181
	idx := 0
182
	for i := 0; i < r; i++ {
183
		for j := 0; j < c; j++ {
184
			if mT.At(i, j) != m.At(i, j) {
185
				return nil, errors.New("matrix is not symmetric")
186
			}
187
			vals[idx] = m.At(i, j)
188
			idx++
189
		}
190
	}
191
	return mat.NewSymDense(r, vals), nil
192
}
193
func RealyPrint(matr *mat.Dense, name string) {
194
	fmatr := mat.Formatted(matr, mat.Prefix(string(strings.Repeat(" ", 2+len(name)))), mat.Squeeze())
195
	fmt.Printf(name+" =%.3v\n", fmatr)
196
}
197

198
func RealyPrint2(matr mat.Dense, name string) {
199
	fmatr := mat.Formatted(&matr, mat.Prefix(string(strings.Repeat(" ", 2+len(name)))), mat.Squeeze())
200
	fmt.Printf(name+" =%.3v\n", fmatr)
201
}
202
func RealyPrintMatrix(matr mat.Matrix, name string) {
203
	fmatr := mat.Formatted(matr, mat.Prefix(string(strings.Repeat(" ", 2+len(name)))), mat.Squeeze())
204
	fmt.Printf(name+" =%.3v\n", fmatr)
205
}
206

207
// Получить медианное значение массива
208
func Median(dataVect mat.VecDense) float64 {
209
	dataVect = SortVecDense(dataVect)
210
	var median float64
211
	l := dataVect.Len()
212
	if l == 0 {
213
		return 0
214
	} else if l%2 == 0 {
215
		median = (dataVect.AtVec(l/2-1) + dataVect.AtVec(l/2)) / 2
216
	} else {
217
		median = dataVect.AtVec(l / 2)
218
	}
219
	return median
220
}
221
func Median_floatArr(dataVect []float64) float64 {
222
	//dataVect = sortVecDense(dataVect)
223
	sort.Float64s(dataVect)
224
	var median float64
225
	l := len(dataVect)
226
	if l == 0 {
227
		return 0
228
	} else if l%2 == 0 {
229
		median = (dataVect[l/2-1] + dataVect[l/2]) / 2
230
	} else {
231
		median = dataVect[l/2]
232
	}
233
	return median
234
}
235

236
// Сортировка вектора массива по возрастанию.
237
func SortVecDense(dataVect mat.VecDense) mat.VecDense {
238
	dataVectLength := dataVect.Len()
239
	for i := 1; i < dataVectLength; i++ {
240
		j := i - 1
241
		for j >= 0 && dataVect.AtVec(j) > dataVect.AtVec(j+1) {
242
			vspom := dataVect.AtVec(j)
243
			dataVect.SetVec(j, dataVect.AtVec(j+1))
244
			dataVect.SetVec(j, vspom)
245
			j--
246
		}
247
	}
248
	return dataVect
249
}
250

251
func Vec_in_ArrFloat(a mat.Vector) []float64 {
252
	b := make([]float64, a.Len())
253
	for i := 0; i < a.Len(); i++ {
254
		b[i] = a.AtVec(i)
255
	}
256
	return b
257
}
258

259
// Диагональ матрицы в зависимости от корреляции k // reference MatLab diag(A,n)
260
func Diag_of_Dense(matr mat.Dense, k int) mat.VecDense {
261
	r, c := matr.Dims()
262
	var matr2 mat.Matrix
263
	if k > 0 {
264
		matr2 = matr.Slice(0, r, k, c)
265
	} else if k < 0 {
266
		matr2 = matr.Slice(-k, r, 0, c)
267
	} else {
268
		matr2 = matr.Slice(0, r, 0, c)
269
	}
270

271
	vect := mat.NewVecDense(mat.DenseCopyOf(matr2).DiagView().Diag(), nil)
272

273
	for i := 0; i < vect.Len(); i++ {
274
		vect.SetVec(i, matr2.At(i, i))
275
	}
276
	return *vect
277
}
278

279
// Поэлементно разделить нулевое значение столбца Matrix на Vector на вектор
280
func Vector_DivElemVec(a mat.Matrix, b mat.Vector) mat.VecDense {
281
	var div_vectors mat.VecDense
282
	var div_Dense mat.Dense
283
	div_Dense.CloneFrom(a)
284
	//fmt.Println(div_Dense.Dims())
285
	asd := div_Dense.ColView(0)
286
	//fmt.Println(">", asd.Len(), b.Len())
287
	div_vectors.DivElemVec(asd, b)
288
	return div_vectors
289
}
290

291
// Предпоследнюю сроку скопировать в последнюю
292
func EditLastRow(matr mat.Dense) mat.Dense {
293
	r, c := matr.Dims()
294

295
	for i := 0; i < c; i++ {
296
		matr.Set(r-1, i, matr.At(r-2, i))
297
	}
298
	return matr
299
}
300