MatrixCalculator/VectorSet.h at main · 321Exusiai/MatrixCalculator · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
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
// =========================================================
// VectorSet.h — 向量组级操作 (Layer 3, 应用层)
// ---------------------------------------------------------
// 职责: 线性无关判定、基提取、秩、Gram-Schmidt 正交化
// 单个向量的原子操作 (加减乘除、点积、范数) 见 vector.h
// =========================================================
#pragma once

#include "vector.h"
#include "matrix.h"  // 已自动引入 RREF.h (见 matrix.h 末尾)
#include<iostream>
#include<vector>
#include<cmath>
#include<stdexcept>
#include<algorithm>

template <typename T>
class VectorSet {
public:
    enum class VectorOrientation {
        Column,  // 默认：教材风格
        Row
    };

private:
    std::vector<Vector<T>> vecs;
    size_t dim;
    size_t rank;
    std::vector<size_t> pivotIndices;
    VectorOrientation orientation;

public:
    VectorSet(const std::vector<Vector<T>>& inputVecs,
              VectorOrientation orient = VectorOrientation::Column)
        : vecs(inputVecs), orientation(orient)
    {
        if (vecs.empty())
            throw std::invalid_argument("Vector set cannot be empty");

        dim = vecs[0].size();
        for (const auto& v : vecs) {
            if (v.size() != dim)
                throw std::invalid_argument("Vectors have inconsistent dimensions");
        }

        Matrix<T> A(
            orientation == VectorOrientation::Column ? dim : vecs.size(),
            orientation == VectorOrientation::Column ? vecs.size() : dim
        );

        if (orientation == VectorOrientation::Column) {
            for (size_t c = 0; c < vecs.size(); c++)
                for (size_t r = 0; r < dim; r++)
                    A.at(r, c) = vecs[c][r];
        } else {
            for (size_t r = 0; r < vecs.size(); r++)
                for (size_t c = 0; c < dim; c++)
                    A.at(r, c) = vecs[r][c];
        }

        RREF<T> rref(std::move(A));
        rref.toRREF();

        rank = rref.getRank();
        pivotIndices = (orient == VectorOrientation::Column)
            ? rref.getPivotCols()
            : rref.getPivotRows();
    }

    // 辅助构造函数：支持原始 std::vector<std::vector<T>>
    VectorSet(const std::vector<std::vector<T>>& inputRaw,
              VectorOrientation orient = VectorOrientation::Column)
        : orientation(orient) {
        for (const auto& row : inputRaw) {
            vecs.emplace_back(row);
        }
        // 不能在构造函数体内赋值委派，直接进行初始化逻辑
        if (vecs.empty())
            throw std::invalid_argument("Vector set cannot be empty");

        dim = vecs[0].size();
        for (const auto& v : vecs) {
            if (v.size() != dim)
                throw std::invalid_argument("Vectors have inconsistent dimensions");
        }

        Matrix<T> A(
            orientation == VectorOrientation::Column ? dim : vecs.size(),
            orientation == VectorOrientation::Column ? vecs.size() : dim
        );

        if (orientation == VectorOrientation::Column) {
            for (size_t c = 0; c < vecs.size(); c++)
                for (size_t r = 0; r < dim; r++)
                    A.at(r, c) = vecs[c][r];
        } else {
            for (size_t r = 0; r < vecs.size(); r++)
                for (size_t c = 0; c < dim; c++)
                    A.at(r, c) = vecs[r][c];
        }

        RREF<T> rref(std::move(A));
        rref.toRREF();

        rank = rref.getRank();
        pivotIndices = (orient == VectorOrientation::Column)
            ? rref.getPivotCols()
            : rref.getPivotRows();
    }

    bool isLinearIndependent() const noexcept {
        return rank == vecs.size();
    }

    std::vector<Vector<T>> basis() const {
        std::vector<Vector<T>> b;
        for (size_t idx : pivotIndices)
            b.push_back(vecs[idx]);
        return b;
    }

    size_t dimension() const noexcept { return rank; }

    static std::vector<Vector<T>>
    gramSchmidt(const std::vector<Vector<T>>& vectors, bool normalize = false) {
        if (vectors.empty())
            throw std::invalid_argument("Gram-Schmidt: empty vector set");

        int dim = vectors[0].size();
        for (const auto& v : vectors)
            if (v.size() != dim)
                throw std::invalid_argument("Gram-Schmidt: dimension mismatch");

        std::vector<Vector<T>> orth;

        for (const auto& v : vectors) {
            Vector<T> u = v;

            for (const auto& uj : orth) {
                T ip_uj = uj.dot(uj);
                if (std::abs(ip_uj) > 1e-9) {
                    T coeff = v.dot(uj) / ip_uj;
                    u -= (uj * coeff);
                }
            }

            if (u.norm() < 1e-9)
                continue;

            if (normalize) {
                u = u.normalized();
            }

            orth.push_back(u);
        }

        return orth;
    }

};