NonlinearFitClass< M, ModelType > Class Template Reference

Collaboration diagram for NonlinearFitClass< M, ModelType >:

Public Member Functions
void	Fit (const Array< double, 1 > &x, const Array< double, 1 > &y, const Array< double, 1 > &sigma, TinyVector< double, M > &params)
Array< double, 2 > &	GetCovariance ()
	NonlinearFitClass (ModelType model)

Private Member Functions
double	Chi2 (const Array< double, 1 > &x, const Array< double, 1 > &y, const Array< double, 1 > &sigma, TinyVector< double, M > params)
void	CalcAlphaBeta (const Array< double, 1 > &x, const Array< double, 1 > &y, const Array< double, 1 > &sigma, TinyVector< double, M > params)
void	Solve ()

Private Attributes
ModelType	Model
TinyVector< double, M >	Params
TinyVector< double, M >	Beta
TinyVector< double, M >	dParams
TinyMatrix< double, M, M >	Alpha
Array< double, 2 >	AlphaInv

Detailed Description

template<int M, typename ModelType>
class NonlinearFitClass< M, ModelType >

Definition at line 24 of file NonlinearFit.h.

Constructor & Destructor Documentation

NonlinearFitClass()

NonlinearFitClass ( ModelType  model )

inline

Definition at line 49 of file NonlinearFit.h.

References NonlinearFitClass< M, ModelType >::AlphaInv, and Array< T, D, ALLOC >::resize().

: Model(model) { AlphaInv.resize(M, M); };

Member Function Documentation

CalcAlphaBeta()

void CalcAlphaBeta ( const Array< double, 1 > &  x, const Array< double, 1 > &  y, const Array< double, 1 > &  sigma, TinyVector< double, M >  params  )

private

Definition at line 77 of file NonlinearFit.h.

References qmcplusplus::Units::force::N, and Array< T, D, ALLOC >::size().

{
  Model.SetParams(params);
  int N = x.size();
  assert(y.size() == N);
  assert(sigma.size() == N);

  for (int k = 0; k < M; k++)
  {
    Beta[k] = 0.0;
    for (int l = 0; l < M; l++)
      Alpha(k, l) = 0.0;
  }

  for (int i = 0; i < N; i++)
  {
    double val                 = Model(x(i));
    TinyVector<double, M> grad = Model.Grad(x(i));
    for (int k = 0; k < M; k++)
      Beta[k] += grad[k] * (y(i) - val) / (sigma(i) * sigma(i));
  }
  for (int i = 0; i < N; i++)
  {
    double val                 = Model(x(i));
    TinyVector<double, M> grad = Model.Grad(x(i));
    for (int k = 0; k < M; k++)
      for (int l = 0; l < M; l++)
        Alpha(k, l) += grad[k] * grad[l] / (sigma(i) * sigma(i));
  }
  // std::cerr << "Alpha = " << Alpha << std::endl;
}

Chi2()

double Chi2 ( const Array< double, 1 > &  x, const Array< double, 1 > &  y, const Array< double, 1 > &  sigma, TinyVector< double, M >  params  )

private

Definition at line 54 of file NonlinearFit.h.

References qmcplusplus::Units::force::N, and Array< T, D, ALLOC >::size().

{
  int N = x.size();
  assert(y.size() == N);
  assert(sigma.size() == N);
  Model.SetParams(params);

  double chi2 = 0;
  for (int i = 0; i < N; i++)
  {
    double val = Model(x(i));
    //     std::cerr << "val = " << val << std::endl;
    //     std::cerr << "yi  = " << y(i) << std::endl;
    //     std::cerr << "sigma = " << sigma(i) << std::endl;
    chi2 += (val - y(i)) * (val - y(i)) / (sigma(i) * sigma(i));
  }
  return chi2;
}

Fit()

void Fit	(	const Array< double, 1 > &	x,
		const Array< double, 1 > &	y,
		const Array< double, 1 > &	sigma,
		TinyVector< double, M > &	params
	)

Definition at line 153 of file NonlinearFit.h.

References BLAS::done, and GJInverse().

{
  double chiNow = Chi2(x, y, sigma, params);
  double lambda = 0.01;

  bool done            = false;
  int iter             = 1;
  int numSmallDecrease = 0;
  while (!done)
  {
    //     std::cerr << "Iteration " << iter << ":  Chi2 = " << chiNow << std::endl;
    //    std::cerr << "params = " << params << std::endl;
    //     std::cerr << "lambda = " << lambda << std::endl;

    CalcAlphaBeta(x, y, sigma, params);
    for (int i = 0; i < M; i++)
      Alpha(i, i) *= (1.0 + lambda);
    Solve();
    //     done = true;
    //     for (int i=0; i<M; i++)
    //       if (std::abs(Beta[i]) > 1.0e-6)
    //  done = false;
    if (!done)
    {
      double chiNew = Chi2(x, y, sigma, params + dParams);
      if (chiNew <= chiNow)
      {
        lambda *= 0.5;
        params = params + dParams;
        if ((chiNow - chiNew) < 0.001)
          numSmallDecrease++;
        if (numSmallDecrease > 4)
          done = true;
        chiNow = chiNew;
      }
      else
      {
        lambda *= 2.0;
      }
    }
    iter++;
  }
  CalcAlphaBeta(x, y, sigma, params);
  for (int i = 0; i < M; i++)
    for (int j = 0; j < M; j++)
      AlphaInv(i, j) = Alpha(i, j);
  GJInverse(AlphaInv);
  //   std::cerr << "Covariace matrix:\n";
  //   for (int i=0; i<AlphaInv.rows(); i++) {
  //     for (int j=0; j<AlphaInv.cols(); j++)
  //       fprintf (stderr, "%12.4e ", AlphaInv(i,j));
  //     fprintf (stderr, "\n");
  //   }
  Model.SetParams(params);
  //  std::cerr << "Chi2 = " << chiNow << std::endl;
}

GetCovariance()

Array<double, 2>& GetCovariance ( )

inline

Definition at line 47 of file NonlinearFit.h.

References NonlinearFitClass< M, ModelType >::AlphaInv.

{ return AlphaInv; }

Solve()

void Solve ( )

private

Definition at line 114 of file NonlinearFit.h.

References GJInverse().

{
  for (int i = 0; i < M; i++)
  {
    dParams[i] = 0.0;
    for (int j = 0; j < M; j++)
      AlphaInv(i, j) = Alpha(i, j);
  }

  GJInverse(AlphaInv);
  TinyMatrix<double, M, M> id;
  id = 0.0;
  for (int i = 0; i < M; i++)
    for (int j = 0; j < M; j++)
      for (int k = 0; k < M; k++)
        id(i, j) += AlphaInv(i, k) * Alpha(k, j);

  for (int i = 0; i < M; i++)
    for (int j = 0; j < M; j++)
      dParams[i] += AlphaInv(i, j) * Beta(j);

  // Do simple Gaussian elimination
  //   double dInv = 1.0/Alpha(0,0);
  //   for (int col=0; col<M; col++)
  //     Alpha(0,col) *= dInv;
  //   Beta(0) *= dInv;
  //   for (int row=1; row<M; row++) {
  //     for (int col=0; col<row; col++) {
  //       Beta(row) -= Alpha(row,col)*Beta(col);
  //       for (int k=0;
  //     double diag = Alpha(row,row);
  //     double diagInv = 1.0/diag;
  //     for (int col=0; col<M; col++)
  //     Beta(row) *= diagInv;
  //     for (int col=row+1; col<
}

Member Data Documentation

Alpha

TinyMatrix<double, M, M> Alpha

private

Definition at line 29 of file NonlinearFit.h.

AlphaInv

Array<double, 2> AlphaInv

private

Definition at line 30 of file NonlinearFit.h.

Referenced by NonlinearFitClass< M, ModelType >::GetCovariance(), and NonlinearFitClass< M, ModelType >::NonlinearFitClass().

Beta

TinyVector<double, M> Beta

private

Definition at line 28 of file NonlinearFit.h.

dParams

TinyVector<double, M> dParams

private

Definition at line 28 of file NonlinearFit.h.

Model

ModelType Model

private

Definition at line 27 of file NonlinearFit.h.

Params

TinyVector<double, M> Params

private

Definition at line 28 of file NonlinearFit.h.

---

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

• /home/pk7/projects/qmc/for_cron_doxygen/qmcpack/src/QMCTools/ppconvert/src/common/NonlinearFit.h