source: ntrip/branches/BNC_2.12/newmat/nm_misc.cpp@ 7868

Last change on this file since 7868 was 2013, checked in by mervart, 13 years ago

* empty log message *

File size: 3.6 KB
Line 
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
14namespace 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
24ReturnMatrix 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
46ReturnMatrix 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
61ReturnMatrix 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
76ReturnMatrix 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
93Real 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
109ReturnMatrix 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
125ReturnMatrix 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///@}
Note: See TracBrowser for help on using the repository browser.