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(n1,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 n1 x n matrix to give n1 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(n1);


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(n1);


69  Y = 0.0;


70  if (j > 1) Y(j1) = sqrt((Real)(j1) / (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 = n1; 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 jth 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 = n1; 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(m1, 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  ///@}

