TSvd Class Reference

#include <xmath.h>

Collaboration diagram for TSvd:

Public Member Functions
void	GetXV (const int RecN, TFltV &VarV) const
double	GetY (const int RecN) const
double	GetSig (const int RecN) const
void	NR_svbksb (TFltVV &u, TFltV &w, TFltVV &v, int m, int n, TFltV &b, TFltV &x)
void	NR_svdvar (TFltVV &v, int ma, TFltV &w, TFltVV &cvm)
void	NR_svdfit ()
	TSvd ()
	~TSvd ()
	TSvd (TSIn &)
void	Save (TSOut &)
TSvd &	operator= (const TSvd &)
int	GetRecs () const
int	GetVars () const
double	GetCf (const int &VarN) const
double	GetCfUncer (const int &VarN) const
double	GetCovar (const int &VarN1, const int &VarN2) const
double	GetChiSq () const
void	GetCfV (TFltV &_CfV)
void	GetCfUncerV (TFltV &CfUncerV)
void	Wr () const

Static Public Member Functions
static double	NR_SIGN (double a, double b)
static double	NR_FMAX (double maxarg1, double maxarg2)
static int	NR_IMIN (int iminarg1, int iminarg2)
static double	NR_pythag (double a, double b)
static void	NR_svdcmp (TFltVV &a, int m, int n, TFltV &w, TFltVV &v)
static PSvd	New (const TFltVV &XVV, const TFltV &YV, const TFltV &SigV=TFltV())
static PSvd	Load (TSIn &SIn)
static void	Svd (const TFltVV &InMtx, TFltVV &LSingV, TFltV &SingValV, TFltVV &RSingV)
static void	Svd1Based (const TFltVV &InMtx1, TFltVV &LSingV, TFltV &SingValV, TFltVV &RSingV)

Public Attributes
TFltVV	XVV
TFltV	YV
TFltV	SigV
int	Recs
int	Vars
TFltVV	CovarVV
TFltV	CfV
double	ChiSq

Private Attributes
TCRef	CRef

Friends
class	TPt< TSvd >

Detailed Description

Definition at line 397 of file xmath.h.

Constructor & Destructor Documentation

TSvd::TSvd ( )

inline

Definition at line 424 of file xmath.h.

Referenced by New().

{}

Here is the caller graph for this function:

TSvd::~TSvd ( )

inline

Definition at line 427 of file xmath.h.

{}

TSvd::TSvd ( TSIn &  )

inline

Definition at line 428 of file xmath.h.

References Fail.

{Fail;}

Member Function Documentation

double TSvd::GetCf ( const int &  VarN ) const

inline

Definition at line 437 of file xmath.h.

Referenced by Wr().

{return CfV[VarN+1];}

Here is the caller graph for this function:

double TSvd::GetCfUncer ( const int &  VarN ) const

inline

Definition at line 438 of file xmath.h.

Referenced by GetCfUncerV(), and Wr().

                                           {
    return sqrt(double(CovarVV.At(VarN+1, VarN+1)));}

Here is the caller graph for this function:

void TSvd::GetCfUncerV

(

TFltV &

CfUncerV

)

Definition at line 1224 of file xmath.cpp.

References TVec< TVal, TSizeTy >::Gen(), GetCfUncer(), and Vars.

                                     {
  CfUncerV.Gen(Vars);
  for (int VarN=0; VarN<Vars; VarN++){
    CfUncerV[VarN]=GetCfUncer(VarN);
  }
}

Here is the call graph for this function:

void TSvd::GetCfV

(

TFltV &

_CfV

)

Definition at line 1220 of file xmath.cpp.

References CfV, and TVec< TVal, TSizeTy >::Del().

                            {
  _CfV=CfV; _CfV.Del(0);
}

Here is the call graph for this function:

double TSvd::GetChiSq ( ) const

inline

Definition at line 443 of file xmath.h.

Referenced by Wr().

{return ChiSq;}

Here is the caller graph for this function:

double TSvd::GetCovar ( const int &  VarN1, const int &  VarN2  ) const

inline

Definition at line 440 of file xmath.h.

Referenced by Wr().

                                                            {
    return CovarVV.At(VarN1+1, VarN2+1);}

Here is the caller graph for this function:

int TSvd::GetRecs ( ) const

inline

Definition at line 434 of file xmath.h.

{return Recs;}

double TSvd::GetSig ( const int  RecN ) const

inline

Definition at line 411 of file xmath.h.

Referenced by NR_svdfit().

{return SigV[RecN-1];}

Here is the caller graph for this function:

int TSvd::GetVars ( ) const

inline

Definition at line 435 of file xmath.h.

{return Vars;}

void TSvd::GetXV ( const int  RecN, TFltV &  VarV  ) const

inline

Definition at line 406 of file xmath.h.

References TVVec< TVal, TSizeTy >::At(), and TVec< TVal, TSizeTy >::Gen().

Referenced by NR_svdfit().

                                                {
    VarV.Gen(Vars+1);
    for (int VarN=0; VarN<Vars; VarN++){VarV[VarN+1]=XVV.At(RecN-1, VarN);}
  }

Here is the call graph for this function:

Here is the caller graph for this function:

double TSvd::GetY ( const int  RecN ) const

inline

Definition at line 410 of file xmath.h.

Referenced by NR_svdfit().

{return YV[RecN-1];}

Here is the caller graph for this function:

static PSvd TSvd::Load ( TSIn &  SIn )

inlinestatic

Definition at line 429 of file xmath.h.

{return new TSvd(SIn);}

PSvd TSvd::New ( const TFltVV &  XVV, const TFltV &  YV, const TFltV &  SigV = TFltV()  )

static

Definition at line 937 of file xmath.cpp.

References TVec< TVal, TSizeTy >::Empty(), IAssert, Svd(), and TSvd().

                                                                      {
  PSvd Svd=PSvd(new TSvd());
  Svd->XVV=_XVV;
  Svd->YV=_YV;
  if (_SigV.Empty()){
    Svd->SigV.Gen(Svd->YV.Len());
    Svd->SigV.PutAll(1);
  } else {
    Svd->SigV=_SigV;
  }
  Svd->Recs=Svd->XVV.GetXDim();
  Svd->Vars=Svd->XVV.GetYDim();
  IAssert(Svd->Recs>0);
  IAssert(Svd->Vars>0);
  IAssert(Svd->YV.Len()==Svd->Recs);
  IAssert(Svd->SigV.Len()==Svd->Recs);
  Svd->CovarVV.Gen(Svd->Vars+1, Svd->Vars+1);
  Svd->CfV.Gen(Svd->Vars+1);
  Svd->NR_svdfit();
  return Svd;
}

Here is the call graph for this function:

static double TSvd::NR_FMAX ( double  maxarg1, double  maxarg2  )

inlinestatic

Definition at line 413 of file xmath.h.

Referenced by NR_svdcmp().

                                                       {
    return maxarg1 > maxarg2 ? maxarg1 : maxarg2;}

Here is the caller graph for this function:

static int TSvd::NR_IMIN ( int  iminarg1, int  iminarg2  )

inlinestatic

Definition at line 415 of file xmath.h.

Referenced by NR_svdcmp().

                                                {
    return iminarg1 < iminarg2 ? iminarg1 : iminarg2;}

Here is the caller graph for this function:

double TSvd::NR_pythag ( double  a, double  b  )

static

Definition at line 959 of file xmath.cpp.

References TMath::Sqr().

Referenced by NR_svdcmp().

                                        {
  double absa,absb;
  absa=fabs(a);
  absb=fabs(b);
  if (absa > absb){
    return absa*sqrt(1.0+TMath::Sqr(absb/absa));
  } else {
    return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+TMath::Sqr(absa/absb)));
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

static double TSvd::NR_SIGN ( double  a, double  b  )

inlinestatic

Definition at line 412 of file xmath.h.

Referenced by NR_svdcmp().

{return b >= 0.0 ? fabs(a) : -fabs(a);}

Here is the caller graph for this function:

void TSvd::NR_svbksb	(	TFltVV &	u,
		TFltV &	w,
		TFltVV &	v,
		int	m,
		int	n,
		TFltV &	b,
		TFltV &	x
	)

Definition at line 1147 of file xmath.cpp.

References TVVec< TVal, TSizeTy >::At().

Referenced by NR_svdfit().

                                                                  {
  int jj,j,i;
  double s;

  TFltV tmp(n+1);
  for (j=1;j<=n;j++) {
    s=0.0;
    if (w[j]) {
      for (i=1;i<=m;i++) s += u.At(i,j)*b[i];
      s /= w[j];
    }
    tmp[j]=s;
  }
  for (j=1;j<=n;j++) {
    s=0.0;
    for (jj=1;jj<=n;jj++) s += v.At(j,jj)*tmp[jj];
    x[j]=s;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void TSvd::NR_svdcmp ( TFltVV &  a, int  m, int  n, TFltV &  w, TFltVV &  v  )

static

Definition at line 970 of file xmath.cpp.

References TVVec< TVal, TSizeTy >::At(), NR_FMAX(), NR_IMIN(), NR_pythag(), NR_SIGN(), and TExcept::Throw().

Referenced by NR_svdfit(), Svd(), and Svd1Based().

                                                                {
  int flag,i,its,j,jj,k,l=0,nm;
  double anorm,c,f,g,h,s,scale,x,y,z;

  TFltV rv1(n+1);
  g=scale=anorm=0.0;
  for (i=1;i<=n;i++) {
    l=i+1;
    rv1[i]=scale*g;
    g=s=scale=0.0;
    if (i <= m) {
      for (k=i;k<=m;k++) scale += fabs(double(a.At(k,i)));
      if (scale) {
        for (k=i;k<=m;k++) {
          a.At(k,i) /= scale;
          s += a.At(k,i)*a.At(k,i);
        }
        f=a.At(i,i);
        g = -NR_SIGN(sqrt(s),f);
        h=f*g-s;
        a.At(i,i)=f-g;
        for (j=l;j<=n;j++) {
          for (s=0.0,k=i;k<=m;k++) s += a.At(k,i)*a(k,j);
          f=s/h;
          for (k=i;k<=m;k++) a.At(k,j) += f*a.At(k,i);
        }
        for (k=i;k<=m;k++) a.At(k,i) *= scale;
      }
    }
    w[i]=scale *g;
    g=s=scale=0.0;
    if (i <= m && i != n) {
      for (k=l;k<=n;k++) scale += fabs(double(a.At(i,k)));
      if (scale) {
        for (k=l;k<=n;k++) {
          a.At(i,k) /= scale;
          s += a.At(i,k)*a.At(i,k);
        }
        f=a.At(i,l);
        g = -NR_SIGN(sqrt(s),f);
        h=f*g-s;
        a.At(i,l)=f-g;
        for (k=l;k<=n;k++) rv1[k]=a.At(i,k)/h;
        for (j=l;j<=m;j++) {
          for (s=0.0,k=l;k<=n;k++) s += a.At(j,k)*a.At(i,k);
          for (k=l;k<=n;k++) a.At(j,k) += s*rv1[k];
        }
        for (k=l;k<=n;k++) a.At(i,k) *= scale;
      }
    }
    anorm=NR_FMAX(anorm,(fabs(double(w[i]))+fabs(double(rv1[i]))));
  }
  for (i=n;i>=1;i--) {
    if (i < n) {
      if (g) {
        for (j=l;j<=n;j++)
          v.At(j,i)=(a.At(i,j)/a.At(i,l))/g;
        for (j=l;j<=n;j++) {
          for (s=0.0,k=l;k<=n;k++) s += a.At(i,k)*v.At(k,j);
          for (k=l;k<=n;k++) v.At(k,j) += s*v.At(k,i);
        }
      }
      for (j=l;j<=n;j++) v.At(i,j)=v.At(j,i)=0.0;
    }
    v.At(i,i)=1.0;
    g=rv1[i];
    l=i;
  }
  for (i=NR_IMIN(m,n);i>=1;i--) {
    l=i+1;
    g=w[i];
    for (j=l;j<=n;j++) a.At(i,j)=0.0;
    if (g) {
      g=1.0/g;
      for (j=l;j<=n;j++) {
        for (s=0.0,k=l;k<=m;k++) s += a.At(k,i)*a.At(k,j);
        f=(s/a.At(i,i))*g;
        for (k=i;k<=m;k++) a.At(k,j) += f*a.At(k,i);
      }
      for (j=i;j<=m;j++) a.At(j,i) *= g;
    } else for (j=i;j<=m;j++) a.At(j,i)=0.0;
    a.At(i,i)++;
  }
  for (k=n;k>=1;k--) {
    for (its=1;its<=30;its++) {
      flag=1;
      for (l=k;l>=1;l--) {
        nm=l-1;
        if ((double)(fabs(double(rv1[l])+anorm)) == anorm) {
          flag=0;
          break;
        }
        if ((double)(fabs(double(w[nm]))+anorm) == anorm) break;
      }
      if (flag) {
        c=0.0;
        s=1.0;
        for (i=l;i<=k;i++) {
          f=s*rv1[i];
          rv1[i]=c*rv1[i];
          if ((double)(fabs(f)+anorm) == anorm) break;
          g=w[i];
          h=NR_pythag(f,g);
          w[i]=h;
          h=1.0/h;
          c=g*h;
          s = -f*h;
          for (j=1;j<=m;j++) {
            y=a.At(j,nm);
            z=a.At(j,i);
            a.At(j,nm)=y*c+z*s;
            a.At(j,i)=z*c-y*s;
          }
        }
      }
      z=w[k];
      if (l == k) {
        if (z < 0.0) {
          w[k] = -z;
          for (j=1;j<=n;j++) v.At(j,k) = -v.At(j,k);
        }
        break;
      }
      if (its==30){
        TExcept::Throw("no convergence in 30 svdcmp iterations");}
      x=w[l];
      nm=k-1;
      y=w[nm];
      g=rv1[nm];
      h=rv1[k];
      f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
      g=NR_pythag(f,1.0);
      f=((x-z)*(x+z)+h*((y/(f+NR_SIGN(g,f)))-h))/x;
      c=s=1.0;
      for (j=l;j<=nm;j++) {
        i=j+1;
        g=rv1[i];
        y=w[i];
        h=s*g;
        g=c*g;
        z=NR_pythag(f,h);
        rv1[j]=z;
        c=f/z;
        s=h/z;
        f=x*c+g*s;
        g = g*c-x*s;
        h=y*s;
        y *= c;
        for (jj=1;jj<=n;jj++) {
          x=v.At(jj,j);
          z=v.At(jj,i);
          v.At(jj,j)=x*c+z*s;
          v.At(jj,i)=z*c-x*s;
        }
        z=NR_pythag(f,h);
        w[j]=z;
        if (z) {
          z=1.0/z;
          c=f*z;
          s=h*z;
        }
        f=c*g+s*y;
        x=c*y-s*g;
        for (jj=1;jj<=m;jj++) {
          y=a.At(jj,j);
          z=a.At(jj,i);
          a.At(jj,j)=y*c+z*s;
          a.At(jj,i)=z*c-y*s;
        }
      }
      rv1[l]=0.0;
      rv1[k]=f;
      w[k]=x;
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void TSvd::NR_svdfit

(

)

Definition at line 1185 of file xmath.cpp.

References TVVec< TVal, TSizeTy >::At(), CfV, ChiSq, CovarVV, TVVec< TVal, TSizeTy >::Gen(), GetSig(), GetXV(), GetY(), NR_svbksb(), NR_svdcmp(), NR_svdvar(), Recs, and Vars.

                    {
  int j,i;
  double wmax,tmp,thresh,sum;
  double TOL=1.0e-5;

  TFltVV u(Recs+1, Vars+1);
  TFltVV v(Vars+1, Vars+1);
  TFltV w(Vars+1);
  TFltV b(Recs+1);
  TFltV afunc(Vars+1);
  for (i=1;i<=Recs;i++) {
    GetXV(i, afunc); // (*funcs)(x[i],afunc,Vars);
    tmp=1.0/GetSig(i);
    for (j=1;j<=Vars;j++){u.At(i,j)=afunc[j]*tmp;}
    b[i]=GetY(i)*tmp;
  }
  NR_svdcmp(u,Recs,Vars,w,v);
  wmax=0.0;
  for (j=1;j<=Vars;j++){
    if (w[j] > wmax){wmax=w[j];}}
  thresh=TOL*wmax;
  for (j=1;j<=Vars;j++){
    if (double(w[j])<thresh){w[j]=0.0;}}
  NR_svbksb(u,w,v,Recs,Vars,b,CfV);
  ChiSq=0.0;
  for (i=1;i<=Recs;i++) {
    GetXV(i, afunc); // (*funcs)(x[i],afunc,Vars);
    for (sum=0.0,j=1;j<=Vars;j++){sum += CfV[j]*afunc[j];}
    ChiSq += (tmp=(GetY(i)-sum)/GetSig(i),tmp*tmp);
  }
  // covariance matrix calculation
  CovarVV.Gen(Vars+1, Vars+1);
  NR_svdvar(v, Vars, w, CovarVV);
}

Here is the call graph for this function:

void TSvd::NR_svdvar	(	TFltVV &	v,
		int	ma,
		TFltV &	w,
		TFltVV &	cvm
	)

Definition at line 1168 of file xmath.cpp.

References TVVec< TVal, TSizeTy >::At().

Referenced by NR_svdfit().

                                                            {
  int k,j,i;
  double sum;

  TFltV wti(ma+1);
  for (i=1;i<=ma;i++) {
    wti[i]=0.0;
    if (w[i]) wti[i]=1.0/(w[i]*w[i]);
  }
  for (i=1;i<=ma;i++) {
    for (j=1;j<=i;j++) {
      for (sum=0.0,k=1;k<=ma;k++) sum += v.At(i,k)*v.At(j,k)*wti[k];
      cvm.At(j,i)=cvm.At(i,j)=sum;
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

TSvd& TSvd::operator= ( const TSvd &  )

inline

Definition at line 432 of file xmath.h.

References Fail.

{Fail; return *this;}

void TSvd::Save ( TSOut &  )

inline

Definition at line 430 of file xmath.h.

References Fail.

{Fail;}

void TSvd::Svd ( const TFltVV &  InMtx, TFltVV &  LSingV, TFltV &  SingValV, TFltVV &  RSingV  )

static

Definition at line 1232 of file xmath.cpp.

References TVVec< TVal, TSizeTy >::At(), TVec< TVal, TSizeTy >::Del(), TVVec< TVal, TSizeTy >::DelX(), TVVec< TVal, TSizeTy >::DelY(), TVec< TVal, TSizeTy >::Gen(), TVVec< TVal, TSizeTy >::Gen(), TVVec< TVal, TSizeTy >::GetXDim(), TVVec< TVal, TSizeTy >::GetYDim(), and NR_svdcmp().

Referenced by TLinAlg::InverseSVD(), New(), and TGraphKey::TakeSig().

                                                                                   {
  //LSingV = InMtx;
  LSingV.Gen(InMtx.GetYDim()+1, InMtx.GetYDim()+1);
  // create 1 based adjacency matrix
  for (int x = 0; x < InMtx.GetXDim(); x++) {
    for (int y = 0; y < InMtx.GetYDim(); y++) {
      LSingV.At(x+1, y+1) = InMtx.At(x, y);
    }
  }
  RSingV.Gen(InMtx.GetYDim()+1, InMtx.GetYDim()+1);
  SingValV.Gen(InMtx.GetYDim()+1);
  TSvd::NR_svdcmp(LSingV, InMtx.GetXDim(), InMtx.GetYDim(), SingValV, RSingV);
  // 0-th singular value/vector is full of zeros, delete it
  SingValV.Del(0);
  LSingV.DelX(0); LSingV.DelY(0);
  RSingV.DelX(0); RSingV.DelY(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

void TSvd::Svd1Based ( const TFltVV &  InMtx1, TFltVV &  LSingV, TFltV &  SingValV, TFltVV &  RSingV  )

static

Definition at line 1252 of file xmath.cpp.

References TVec< TVal, TSizeTy >::Del(), TVVec< TVal, TSizeTy >::DelX(), TVVec< TVal, TSizeTy >::DelY(), TVec< TVal, TSizeTy >::Gen(), TVVec< TVal, TSizeTy >::Gen(), TVVec< TVal, TSizeTy >::GetXDim(), TVVec< TVal, TSizeTy >::GetYDim(), and NR_svdcmp().

Referenced by TSnap::GetSngVals(), and TSnap::GetSngVec().

                                                                                          {
  LSingV = InMtx1;
  SingValV.Gen(InMtx1.GetYDim());
  RSingV.Gen(InMtx1.GetYDim(), InMtx1.GetYDim());
  TSvd::NR_svdcmp(LSingV, InMtx1.GetXDim()-1, InMtx1.GetYDim()-1, SingValV, RSingV);
  // 0-th singular value/vector is full of zeros, delete it
  SingValV.Del(0);
  LSingV.DelX(0); LSingV.DelY(0);
  RSingV.DelX(0); RSingV.DelY(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

void TSvd::Wr

(

)

const

Definition at line 1263 of file xmath.cpp.

References GetCf(), GetCfUncer(), GetChiSq(), GetCovar(), and Vars.

                    {
  printf("\n%11s %21s\n","parameter","uncertainty");
  for (int i=0; i<Vars;i++){
    printf("  a[%1d] = %8.6f %12.6f\n",
     i+1, GetCf(i), GetCfUncer(i));
  }
  printf("chi-squared = %12f\n", GetChiSq());

  printf("full covariance matrix\n");
  {for (int i=0;i<Vars;i++){
    for (int j=0;j<Vars;j++){
      printf("%12f", GetCovar(i, j));}
    printf("\n");
  }}
}

Here is the call graph for this function:

Friends And Related Function Documentation

friend class TPt< TSvd >

friend

Definition at line 397 of file xmath.h.

Member Data Documentation

TFltV TSvd::CfV

Definition at line 404 of file xmath.h.

Referenced by GetCfV(), and NR_svdfit().

double TSvd::ChiSq

Definition at line 405 of file xmath.h.

Referenced by NR_svdfit().

TFltVV TSvd::CovarVV

Definition at line 403 of file xmath.h.

Referenced by NR_svdfit().

TCRef TSvd::CRef

private

Definition at line 397 of file xmath.h.

int TSvd::Recs

Definition at line 402 of file xmath.h.

Referenced by NR_svdfit().

TFltV TSvd::SigV

Definition at line 401 of file xmath.h.

int TSvd::Vars

Definition at line 402 of file xmath.h.

Referenced by GetCfUncerV(), NR_svdfit(), and Wr().

TFltVV TSvd::XVV

Definition at line 399 of file xmath.h.

TFltV TSvd::YV

Definition at line 400 of file xmath.h.

---

The documentation for this class was generated from the following files:

• glib-core/xmath.h
• glib-core/xmath.cpp