ssa

1
package main
2

3
import (
4
	"fmt"
5
	"log"
6

7
	"github.com/RB-PRO/ssa/pkg/oss"
8
	"gonum.org/v1/gonum/mat"
9
)
10

11
type SDVs struct {
12
	V *mat.Dense
13
	U *mat.Dense
14
	S *mat.Dense
15
}
16

17
func AaT(matr *mat.Dense) *mat.Dense { // Multipy matrix AT*A
18
	a := mat.Matrix(matr)
19
	aT := a.T()
20
	ad := mat.DenseCopyOf(a)
21
	aTd := mat.DenseCopyOf(aT)
22
	n1, _ := aTd.Dims()
23
	_, m2 := ad.Dims()
24
	output := mat.NewDense(n1, m2, nil)
25
	fmt.Print("X: ")
26
	fmt.Println(ad.Dims())
27
	fmt.Print("XT: ")
28
	fmt.Println(aTd.Dims())
29
	// oss.SafeToXlsxMatrix(ad, "ad")
30
	// oss.SafeToXlsxMatrix(aTd, "aTd")
31
	output.Mul(ad, aTd)
32
	return output
33
}
34
func SDV(X *mat.Dense, rank int) []SDVs {
35
	SDVsout := make([]SDVs, rank)
36
	for i := 0; i < rank; i++ {
37
		/*
38
			// Это с дроблением на матрицы
39
			X_y, _ := X.Dims()
40
			kk := X.Slice(0, X_y, i, i+X_y) // Взять часть матрицы X
41
			kek := mat.DenseCopyOf(kk)      // Преобразовать в Dense
42
			SDVsout[i] = SDV_single(kek)    // Сохранить значение
43
		*/
44

45
		// это без дробления на матрицы
46
		SDVsout[i] = SDV_single(X)
47
	}
48
	return SDVsout
49
}
50
func SDV_single(matT *mat.Dense) SDVs {
51
	// oss.SafeToXlsxMatrix(matT, "matT")
52

53
	var SDVout SDVs
54
	var svdMat mat.SVD
55
	ok := svdMat.Factorize(matT, mat.SVDFull)
56
	if !ok {
57
		log.Fatal("failed to factorize A")
58
	}
59

60
	SDVout.V = new(mat.Dense)
61
	SDVout.S = new(mat.Dense)
62
	svdMat.VTo(SDVout.V)
63

64
	svdMat.UTo(SDVout.S)
65
	lenX_s, lenY_s := matT.Dims()
66
	//fmt.Println(lenY_s)
67
	valuesMat := make([]float64, lenX_s)
68
	//fmt.Println(len(valuesMat))
69
	svdMat.Values(valuesMat)
70

71
	SDVout.U = mat.NewDense(lenX_s, lenY_s, nil)
72
	for ind, val := range valuesMat {
73
		SDVout.U.Set(ind, ind, val)
74
	}
75

76
	//fmt.Println(SDVout.S.Dims())
77

78
	SDVout.U = mat.DenseCopyOf(SDVout.U.T())
79

80
	//fmt.Println(SDVout.S.Dims())
81

82
	return SDVout
83
}
84

85
func makeRank(matr *mat.Dense) int {
86
	var svd mat.SVD
87
	ok := svd.Factorize(matr, mat.SVDFull)
88
	if !ok {
89
		log.Fatal("failed to factorize A")
90
	}
91
	rank := svd.Rank(oss.Rcond)
92
	if rank == 0 {
93
		log.Fatal("zero rank system")
94
	}
95
	return (rank)
96
}
97