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 | ///@}
|
---|