nm_misc.cpp
Go to the documentation of this file.
1 
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 
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 
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 
ReturnMatrix Helmert(int n, bool full)
Definition: nm_misc.cpp:24
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
ReturnMatrix for_return() const
Definition: newmat4.cpp:249
ReturnMatrix Helmert_transpose(const ColumnVector &Y, bool full)
Definition: nm_misc.cpp:76
double Real
Definition: include.h:307
#define REPORT
Definition: nm_misc.cpp:21
void resize(int)
Definition: newmat4.cpp:318
GetSubMatrix Column(int f) const
Definition: newmat.h:2152
#define Throw(E)
Definition: myexcept.h:191
The usual rectangular matrix.
Definition: newmat.h:625
Index exception.
Definition: newmat.h:1960
int ncols() const
Definition: newmat.h:500
virtual void resize(int, int)
Definition: newmat4.cpp:289
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
int nrows() const
Definition: newmat.h:499
Column vector.
Definition: newmat.h:1008
GetSubMatrix submatrix(int, int, int, int) const
Definition: submat.cpp:27
GetSubMatrix row(int) const
Definition: submat.cpp:49
void release()
Definition: newmat.h:517


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16