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
Helmert(int n,bool full)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
Helmert(const ColumnVector & X,bool full)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
Helmert(int n,int j,bool full)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
Helmert_transpose(const ColumnVector & Y,bool full)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
Helmert_transpose(const ColumnVector & Y,int j,bool full)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
Helmert(const Matrix & X,bool full)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
Helmert_transpose(const Matrix & Y,bool full)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 ///@}
148