1 | /// \ingroup newmat |
---|
2 | ///@{ |
---|
3 | |
---|
4 | /// \file nm_misc.cpp |
---|
5 | /// Miscellaneous programs. |
---|
6 | |
---|
7 | #define WANT_MATH |
---|
8 | |
---|
9 | #include "include.h" |
---|
10 | |
---|
11 | #include "newmatap.h" |
---|
12 | |
---|
13 | #ifdef use_namespace |
---|
14 | namespace NEWMAT { |
---|
15 | #endif |
---|
16 | |
---|
17 | |
---|
18 | #ifdef DO_REPORT |
---|
19 | #define REPORT { static ExeCounter ExeCount(__LINE__,21); ++ExeCount; } |
---|
20 | #else |
---|
21 | #define REPORT {} |
---|
22 | #endif |
---|
23 | |
---|
24 | ReturnMatrix Helmert(int n, bool full) |
---|
25 | { |
---|
26 | REPORT |
---|
27 | Tracer et("Helmert "); |
---|
28 | if (n <= 0) Throw(ProgramException("Dimension <= 0 ")); |
---|
29 | Matrix H; |
---|
30 | |
---|
31 | if (full) H.resize(n,n); else H.resize(n-1,n); |
---|
32 | H = 0.0; |
---|
33 | for (int i = 1; i < n; ++i) |
---|
34 | { |
---|
35 | Real f = 1.0 / sqrt((Real)i * (i+1)); |
---|
36 | H.submatrix(i,i,1,i) = -f; H(i,i+1) = f * i; |
---|
37 | } |
---|
38 | if (full) { H.row(n) = 1.0 / sqrt((Real)n); } |
---|
39 | H.release(); return H.for_return(); |
---|
40 | } |
---|
41 | |
---|
42 | |
---|
43 | |
---|
44 | // Multiply X by n-1 x n matrix to give n-1 contrasts |
---|
45 | // Return a ColumnVector |
---|
46 | ReturnMatrix Helmert(const ColumnVector& X, bool full) |
---|
47 | { |
---|
48 | REPORT |
---|
49 | Tracer et("Helmert * CV"); |
---|
50 | int n = X.nrows(); |
---|
51 | if (n == 0) Throw(ProgramException("X Vector of length 0", X)); |
---|
52 | Real sum = 0.0; ColumnVector Y; |
---|
53 | if (full) Y.resize(n); else Y.resize(n-1); |
---|
54 | for (int i = 1; i < n; ++i) |
---|
55 | { sum += X(i); Y(i) = (i * X(i+1) - sum) / sqrt((Real)i * (i+1)); } |
---|
56 | if (full) { sum += X(n); Y(n) = sum / sqrt((Real)n); } |
---|
57 | Y.release(); return Y.for_return(); |
---|
58 | } |
---|
59 | |
---|
60 | // same as above for X a ColumnVector, length n, element j = 1; otherwise 0 |
---|
61 | ReturnMatrix Helmert(int n, int j, bool full) |
---|
62 | { |
---|
63 | REPORT |
---|
64 | Tracer et("Helmert:single element "); |
---|
65 | if (n <= 0) Throw(ProgramException("X Vector of length <= 0")); |
---|
66 | if (j > n || j <= 0) |
---|
67 | Throw(ProgramException("Out of range element number ")); |
---|
68 | ColumnVector Y; if (full) Y.resize(n); else Y.resize(n-1); |
---|
69 | Y = 0.0; |
---|
70 | if (j > 1) Y(j-1) = sqrt((Real)(j-1) / (Real)j); |
---|
71 | for (int i = j; i < n; ++i) Y(i) = - 1.0 / sqrt((Real)i * (i+1)); |
---|
72 | if (full) Y(n) = 1.0 / sqrt((Real)n); |
---|
73 | Y.release(); return Y.for_return(); |
---|
74 | } |
---|
75 | |
---|
76 | ReturnMatrix Helmert_transpose(const ColumnVector& Y, bool full) |
---|
77 | { |
---|
78 | REPORT |
---|
79 | Tracer et("Helmert_transpose * CV "); |
---|
80 | int n = Y.nrows(); Real sum; |
---|
81 | if (full) sum = Y(n) / sqrt((Real)n); else { sum = 0.0; ++n; } |
---|
82 | ColumnVector X(n); |
---|
83 | for (int i = n-1; i > 0; --i) |
---|
84 | { |
---|
85 | Real Yi = Y(i) / sqrt((Real)i * (i+1)); |
---|
86 | X(i+1) = i * Yi + sum; sum -= Yi; |
---|
87 | } |
---|
88 | X(1) = sum; |
---|
89 | X.release(); return X.for_return(); |
---|
90 | } |
---|
91 | |
---|
92 | // same as above but want only j-th element |
---|
93 | Real Helmert_transpose(const ColumnVector& Y, int j, bool full) |
---|
94 | { |
---|
95 | REPORT |
---|
96 | Tracer et("Helmert_transpose:single element "); |
---|
97 | int n = Y.nrows(); Real sum; |
---|
98 | if (full) sum = Y(n) / sqrt((Real)n); else { sum = 0.0; ++n; } |
---|
99 | if (j > n || j <= 0) Throw(IndexException(j, Y)); |
---|
100 | for (int i = n-1; i > 0; --i) |
---|
101 | { |
---|
102 | Real Yi = Y(i) / sqrt((Real)i * (i+1)); |
---|
103 | if (i+1 == j) return i * Yi + sum; |
---|
104 | sum -= Yi; |
---|
105 | } |
---|
106 | return sum; |
---|
107 | } |
---|
108 | |
---|
109 | ReturnMatrix Helmert(const Matrix& X, bool full) |
---|
110 | { |
---|
111 | REPORT |
---|
112 | Tracer et("Helmert * Matrix"); |
---|
113 | int m = X.nrows(); int n = X.ncols(); |
---|
114 | if (m == 0) Throw(ProgramException("Matrix has 0 rows ", X)); |
---|
115 | Matrix Y; |
---|
116 | if (full) Y.resize(m,n); else Y.resize(m-1, n); |
---|
117 | for (int j = 1; j <= n; ++j) |
---|
118 | { |
---|
119 | ColumnVector CV = X.Column(j); |
---|
120 | Y.Column(j) = Helmert(CV, full); |
---|
121 | } |
---|
122 | Y.release(); return Y.for_return(); |
---|
123 | } |
---|
124 | |
---|
125 | ReturnMatrix Helmert_transpose(const Matrix& Y, bool full) |
---|
126 | { |
---|
127 | REPORT |
---|
128 | Tracer et("Helmert_transpose * Matrix "); |
---|
129 | int m = Y.nrows(); int n = Y.ncols(); if (!full) ++m; |
---|
130 | Matrix X(m, n); |
---|
131 | for (int j = 1; j <= n; ++j) |
---|
132 | { |
---|
133 | ColumnVector CV = Y.Column(j); |
---|
134 | X.Column(j) = Helmert_transpose(CV, full); |
---|
135 | } |
---|
136 | X.release(); return X.for_return(); |
---|
137 | } |
---|
138 | |
---|
139 | |
---|
140 | |
---|
141 | |
---|
142 | #ifdef use_namespace |
---|
143 | } |
---|
144 | #endif |
---|
145 | |
---|
146 | |
---|
147 | ///@} |
---|