home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / OLSSVD.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  5.6 KB  |  226 lines
  1. /**************************************************/
  2. /*        olssvd.c source for olsSvd class        */
  3. /**************************************************/
  4.  
  5.  
  6. /**************************************************/
  7. /*            MatClass Source File                */
  8. /*       Copyright of C. R. Birchenhall           */
  9. /*       University of Manchester, UK.            */
  10. /*   MatClass is freeware. This file should be    */
  11. /* made freely available to users of any software */
  12. /* whose creation is wholly or partly dependent   */
  13. /*                on this file.                   */
  14. /**************************************************/
  15.  
  16. #include "olssvd.hpp"
  17. #include "matsvd.hpp"
  18. #include <math.h>
  19.  
  20. /*************************************************************/
  21. /*                      olsSvd methods                       */
  22. /*************************************************************/
  23.  
  24. #ifdef __cplusplus
  25. olsSvd::olsSvd( void ) : matOls( )
  26. #else
  27. olsSvd::olsSvd( void ) : ()
  28. #endif
  29. {
  30.    status = 0 ;
  31.    static char *mName = "olsSvd" ;
  32.    matFunc func( mName ) ; 
  33.    SV.name( "olsSV" ) ;
  34.    U.name( "olsU" ) ;
  35.    W.name( "olsW" ) ;
  36.    func.trace(TRUE) ;
  37. } // olsSvd
  38.  
  39. #ifdef __cplusplus
  40. olsSvd::olsSvd( olsSvd& ols ) : matOls( ols )
  41. #else
  42. olsSvd::olsSvd( olsSvd& ols  ) : ( ols )
  43. #endif
  44. {} // olsSvd
  45.  
  46. void olsSvd::operator = ( olsSvd& ols ) 
  47. {
  48.       matOls::operator = ( ols ) ;
  49. } // olsSvd =
  50.  
  51. #ifdef __cplusplus
  52. olsSvd::olsSvd( const matrix& y, const matrix& x )
  53.      : matOls( y, x )
  54. #else
  55. olsSvd::olsSvd( const matrix& y, const matrix& x )
  56.      : ( y, x )
  57. #endif
  58. {
  59.    static char *mName = "olsSvd(y,x)" ;
  60.    matFunc func( mName ) ; 
  61.    SV.name( "olsSV" ) ;
  62.    U.name( "olsU" ) ;
  63.    W.name( "olsW" ) ;
  64.    func.trace(TRUE) ;
  65. } // olsSvd( matrix& )
  66.  
  67. olsSvd::~olsSvd( void )
  68. {
  69.    static char *mName = "~olsSvd" ;
  70.    matFunc func( mName ) ; debugInfo( func ) ;
  71. } // ~olsSvd
  72.  
  73. void olsSvd::initial( void )
  74. {
  75.    matOls::initial() ;
  76.    SV.reset( nVars ) ;
  77.    U.reset( nObs, nVars ) ;
  78.    W.reset( nVars, nVars ) ;
  79. } // olsSvd::initial
  80.  
  81. void olsSvd::clear( void )
  82. {
  83.    matOls::clear() ;
  84.    SV.clear() ;
  85.    U.clear() ;
  86.    W.clear() ;
  87. } // olsSvd::clear
  88.  
  89. void olsSvd::zeroSV( void )
  90. {
  91.    static char *mName = "zero" ;
  92.    matFunc func( mName ) ; debugInfo( func ) ;
  93.    INDEX  i, n = SV.nRows() ;
  94.    REAL   svLimit ;
  95.    for ( svLimit = SV(1), i = 2 ; i <= n ; i++ ) {
  96.       if ( SV(i) > svLimit )
  97.          svLimit = SV(i) ;
  98.    } // for
  99.    svLimit *= tol ;
  100.    for ( i = SV.nRows() ; i ; i-- ) {
  101.       if ( SV(i) < svLimit )
  102.      SV(i) = 0.0 ;
  103.    } // for
  104. } // olsSvd zero
  105.  
  106. void olsSvd::formBeta( void )
  107. {
  108.    static char *mName = "formBeta" ;
  109.    matFunc func( mName ) ; debugInfo( func ) ;
  110.    matError error ;
  111.    refMatrix yRef( Y ), bRef( beta ) ;
  112.    INDEX nc = Y.nCols(), j ;
  113.    for ( j = 1 ; j <= nc ; j++ ) {
  114.        yRef.refCol(j) ;
  115.        bRef.refCol(j) ;
  116.        svdBackSub( U, SV, W, yRef, bRef, error ) ;
  117.        if ( error )
  118.       errorExit( mName, error ) ;
  119.    } // for 
  120. } // olsSvd formBeta
  121.  
  122. void olsSvd::decompose( void )
  123. {
  124.    static char *mName = "decomp" ;
  125.    matFunc func( mName ) ; debugInfo( func ) ;
  126.    matError error ;
  127.    if ( !( status & ASSIGNED ) )
  128.       errorExit( mName, UNASS ) ;
  129.    if ( !( status & DECOMPOSED ) ) {
  130.       U = X ;
  131.       if ( !svdcmp( U, SV, W, error ) )
  132.          errorExit( mName, error ) ;
  133.       zeroSV() ;
  134.       formBeta() ;
  135.       status |= DECOMPOSED ;
  136.    } // if
  137. } // olsSvd::decompose
  138.  
  139. matrix& olsSvd::sv( matrix& s )
  140. {
  141.    static char *mName = "sv" ;
  142.    matFunc func( mName ) ; debugInfo( func ) ;
  143.    decompose() ;
  144.    s = SV ;
  145.    return s ;
  146. } // olsSvd sv
  147.  
  148. void olsSvd::setSV( const matrix& newSV )
  149. {
  150.    static char *mName = "sv" ;
  151.    matFunc func( mName ) ; debugInfo( func ) ;
  152.    decompose() ;
  153.    SV = newSV ;
  154.    formBeta() ;
  155.    formV() ;
  156. } // olsSvd setSV
  157.  
  158. outFile& olsSvd::info( outFile& f ) M_CONST
  159. {
  160.    return matOls::olsInfo( f, "olsSvd" ) ;
  161. } // olsSvd::info
  162.  
  163. REAL olsSvd::varAdd( const matrix& z, INDEX i )
  164. /***************************************************
  165.   Return modified RSS when vars z added to nth eq.
  166. ***************************************************/
  167. {
  168.    static char *mName = "addVar" ;
  169.    matFunc func( mName ) ; debugInfo( func ) ;
  170.    z.errorExit( mName, NIMPL ) ;
  171.    return (REAL) i ;
  172. } // olsSvd::varAdd
  173.  
  174.  
  175. REAL olsSvd::cond( void )
  176. /************************************
  177.    Condition estimate based on SVD
  178. ************************************/
  179. {
  180.    static char *mName = "cond" ;
  181.    matFunc func( mName ) ; debugInfo( func ) ;
  182.    decompose() ;
  183.    if ( !(status & CONDITIONED ) ) {
  184.       REAL min, max ;
  185.       min = max = SV(1) ;
  186.       for ( INDEX i = 2; i <= nVars ; i++ ) {
  187.          if ( SV(i) < min )
  188.             min = SV(i) ;
  189.          if ( SV(i) > max )
  190.             max = SV(i) ;
  191.       } // for i
  192.       if ( min < tol )
  193.          condition = -1 ;
  194.       else
  195.          condition = max/min ;
  196.    } // if 
  197.    return condition ;
  198. } // olsSvd::cond
  199.  
  200.  
  201. void olsSvd::formV( void )
  202. {
  203.    static char *mName = "formV" ;
  204.    matFunc func( mName ) ; debugInfo( func ) ;
  205.    formResid() ;
  206.    if ( !(status & VMATRIX ) ) {
  207.       INDEX n = nVars, i ;
  208.       REAL r ;
  209.       matrix d(n) ;
  210.       for ( i = 1 ; i <= n ; i++ ) {
  211.           r = SV(i) ;
  212.           if ( r!= 0.0 )
  213.              d(i) = 1.0 / ( r * r ) ;
  214.           else
  215.              d(i) = 0.0 ;
  216.       } // for
  217.       // R = W * diag(d) ;
  218.       R.colMultOf( W, d ) ;
  219.       // V = R * W' = W * diag(d) * W'
  220.       V.multTOf( W, R ) ;
  221.       for ( i = 1 ; i <= n ; i++ )
  222.           VSqrt(i) = sqrt( V(i,i) ) ;
  223.       status |= VMATRIX ;
  224.    } // if
  225. } // olsSvd::formV
  226.