#include <xmath.h>

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 395 of file xmath.h.

Constructor & Destructor Documentation

TSvd::TSvd ( ) [inline]

Definition at line 422 of file xmath.h.

{}

TSvd::~TSvd ( ) [inline]

Definition at line 425 of file xmath.h.

{}

TSvd::TSvd ( TSIn & ) [inline]

Definition at line 426 of file xmath.h.

{Fail;}

Member Function Documentation

double TSvd::GetCf ( const int & VarN ) const [inline]

Definition at line 435 of file xmath.h.

{return CfV[VarN+1];}

double TSvd::GetCfUncer ( const int & VarN ) const [inline]

Definition at line 436 of file xmath.h.

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

void TSvd::GetCfUncerV ( TFltV & CfUncerV )

Definition at line 1213 of file xmath.cpp.

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

void TSvd::GetCfV ( TFltV & _CfV )

Definition at line 1209 of file xmath.cpp.

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

double TSvd::GetChiSq ( ) const [inline]

Definition at line 441 of file xmath.h.

{return ChiSq;}

double TSvd::GetCovar	(	const int &	VarN1,
		const int &	VarN2
	)		const `[inline]`

Definition at line 438 of file xmath.h.

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

int TSvd::GetRecs ( ) const [inline]

Definition at line 432 of file xmath.h.

{return Recs;}

double TSvd::GetSig ( const int RecN ) const [inline]

Definition at line 409 of file xmath.h.

{return SigV[RecN-1];}

int TSvd::GetVars ( ) const [inline]

Definition at line 433 of file xmath.h.

{return Vars;}

void TSvd::GetXV	(	const int	RecN,
		TFltV &	VarV
	)		const `[inline]`

Definition at line 404 of file xmath.h.

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

double TSvd::GetY ( const int RecN ) const [inline]

Definition at line 408 of file xmath.h.

{return YV[RecN-1];}

static PSvd TSvd::Load ( TSIn & SIn ) [inline, static]

Definition at line 427 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 926 of file xmath.cpp.

                                                                      {
  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;
}

static double TSvd::NR_FMAX	(	double	maxarg1,
		double	maxarg2
	)		`[inline, static]`

Definition at line 411 of file xmath.h.

                                                       {
    return maxarg1 > maxarg2 ? maxarg1 : maxarg2;}

static int TSvd::NR_IMIN	(	int	iminarg1,
		int	iminarg2
	)		`[inline, static]`

Definition at line 413 of file xmath.h.

                                                {
    return iminarg1 < iminarg2 ? iminarg1 : iminarg2;}

double TSvd::NR_pythag	(	double	a,
		double	b
	)		`[static]`

Definition at line 948 of file xmath.cpp.

                                        {
  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)));
  }
}

static double TSvd::NR_SIGN	(	double	a,
		double	b
	)		`[inline, static]`

Definition at line 410 of file xmath.h.

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

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

Definition at line 1136 of file xmath.cpp.

                                                                  {
  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;
  }
}

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

Definition at line 959 of file xmath.cpp.

                                                                {
  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;
    }
  }
}

void TSvd::NR_svdfit ( )

Definition at line 1174 of file xmath.cpp.

                    {
  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);
}

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

Definition at line 1157 of file xmath.cpp.

                                                            {
  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;
    }
  }
}

TSvd& TSvd::operator= ( const TSvd & ) [inline]

Definition at line 430 of file xmath.h.

{Fail; return *this;}

void TSvd::Save ( TSOut & ) [inline]

Definition at line 428 of file xmath.h.

{Fail;}

void TSvd::Svd	(	const TFltVV &	InMtx,
		TFltVV &	LSingV,
		TFltV &	SingValV,
		TFltVV &	RSingV
	)		`[static]`

Definition at line 1221 of file xmath.cpp.

                                                                                   {
  //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);
}

void TSvd::Svd1Based	(	const TFltVV &	InMtx1,
		TFltVV &	LSingV,
		TFltV &	SingValV,
		TFltVV &	RSingV
	)		`[static]`

Definition at line 1241 of file xmath.cpp.

                                                                                          {
  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);
}

void TSvd::Wr ( ) const

Definition at line 1252 of file xmath.cpp.

                    {
  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");
  }}
}