JCALIBRATE::JFit Class Reference

Fit. More...

#include <JFitK40.hh>

Classes
struct	result_type
	Result type. More...

Public Types
typedef std::shared_ptr < JMEstimator >	estimator_type

Public Member Functions
	JFit (const int option, const int debug)
	Constructor. More...
result_type	operator() (const data_type &data)
	Fit. More...

Public Attributes
int	debug
estimator_type	estimator
	M-Estimator function. More...
double	lambda
JModel	value
int	numberOfIterations
JMATH::JMatrixNS	V

Static Public Attributes
static constexpr int	MAXIMUM_ITERATIONS = 10000
	maximal number of iterations. More...
static constexpr double	EPSILON = 1.0e-4
	maximal distance to minimum. More...
static constexpr double	LAMBDA_MIN = 0.01
	minimal value control parameter More...
static constexpr double	LAMBDA_MAX = 100.0
	maximal value control parameter More...
static constexpr double	LAMBDA_UP = 10.0
	multiplication factor control parameter More...
static constexpr double	LAMBDA_DOWN = 10.0
	multiplication factor control parameter More...
static constexpr double	PIVOT = std::numeric_limits<double>::epsilon()
	minimal value diagonal element of matrix More...

Private Member Functions
void	evaluate (const data_type &data)
	Evaluation of fit. More...

Private Attributes
JMATH::JVectorND	Y
double	successor
JModel	previous
std::vector< double >	h

Detailed Description

Fit.

Definition at line 1118 of file JFitK40.hh.

Member Typedef Documentation

typedef std::shared_ptr<JMEstimator> JCALIBRATE::JFit::estimator_type

Definition at line 1129 of file JFitK40.hh.

Constructor & Destructor Documentation

JCALIBRATE::JFit::JFit ( const int  option, const int  debug  )

inline

Constructor.

Parameters

option	M-estimator
debug	debug

Definition at line 1138 of file JFitK40.hh.

                                            :
      debug(debug)
    {
      using namespace JPP;

      estimator.reset(getMEstimator(option));
    }

Member Function Documentation

result_type JCALIBRATE::JFit::operator() ( const data_type &  data )

inline

Fit.

Parameters

data

data

Returns

chi2, NDF

Definition at line 1153 of file JFitK40.hh.

    {
      using namespace std;
      using namespace JPP;


      value.setIndex();

      const size_t N = value.getN();

      V.resize(N);
      Y.resize(N);
      h.resize(N);

      int ndf = 0;

      for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {

        const pair_type& pair = ix->first;

        if (value.parameters[pair.first ].status &&
            value.parameters[pair.second].status) {

          ndf += ix->second.size();
        }
      }

      ndf -= value.getN();


      lambda    = LAMBDA_MIN;

      double precessor = numeric_limits<double>::max();

      for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {

        DEBUG("step:     " << numberOfIterations << endl);

        evaluate(data);

        DEBUG("lambda:   " << FIXED(12,5) << lambda    << endl);
        DEBUG("chi2:     " << FIXED(12,3) << successor << endl);

        if (successor < precessor) {

          if (numberOfIterations != 0) {

            if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
              return { successor / estimator->getRho(1.0), ndf };
            }

            if (lambda > LAMBDA_MIN) {
              lambda /= LAMBDA_DOWN;
            }
          }

          precessor = successor;
          previous  = value;

        } else {

          value   = previous;
          lambda *= LAMBDA_UP;

          if (lambda > LAMBDA_MAX) {
            return { precessor / estimator->getRho(1.0), ndf }; // no improvement found
          }

          evaluate(data);
        }

        if (debug >= debug_t) {

          size_t row = 0;

          if (value.R .isFree()) { cout << "R  " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.p1.isFree()) { cout << "p1 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.p2.isFree()) { cout << "p2 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.p3.isFree()) { cout << "p3 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.p4.isFree()) { cout << "p4 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.bg.isFree()) { cout << "bg " << FIXED(12,3) << Y[row] << endl; ++row; }
          if (value.cc.isFree()) { cout << "cc " << FIXED(12,3) << Y[row] << endl; ++row; }

          for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
            if (value.parameters[pmt].QE .isFree()) { cout << "PMT[" << setw(2) << pmt << "].QE  " << FIXED(12,5) << Y[row] << endl; ++row; }
            if (value.parameters[pmt].TTS.isFree()) { cout << "PMT[" << setw(2) << pmt << "].TTS " << FIXED(12,5) << Y[row] << endl; ++row; }
            if (value.parameters[pmt].t0 .isFree()) { cout << "PMT[" << setw(2) << pmt << "].t0  " << FIXED(12,5) << Y[row] << endl; ++row; }
          }
        }

        // force definite positiveness

        for (size_t i = 0; i != N; ++i) {

          if (V(i,i) < PIVOT) {
            V(i,i) = PIVOT;
          }

          h[i] = 1.0 / sqrt(V(i,i));
        }
 
        // normalisation

        for (size_t i = 0; i != N; ++i) {
          for (size_t j = 0; j != i; ++j) {
            V(j,i) *= h[i] * h[j];
            V(i,j)  = V(j,i);
          }
        }

        for (size_t i = 0; i != N; ++i) {
          V(i,i) = 1.0 + lambda;
        }

        // solve A x = b

        for (size_t col = 0; col != N; ++col) {
          Y[col] *= h[col];
        }

        try {
          V.solve(Y);
        }
        catch (const exception& error) {

          ERROR("JGandalf: " << error.what() << endl << V << endl);

          break;
        }

        // update value

        const double factor = 2.0;

        size_t row = 0;

        if (value.R .isFree()) { value.R  -= factor * h[row] * Y[row];  ++row; }
        if (value.p1.isFree()) { value.p1 -= factor * h[row] * Y[row];  ++row; }
        if (value.p2.isFree()) { value.p2 -= factor * h[row] * Y[row];  ++row; }
        if (value.p3.isFree()) { value.p3 -= factor * h[row] * Y[row];  ++row; }
        if (value.p4.isFree()) { value.p4 -= factor * h[row] * Y[row];  ++row; }
        if (value.bg.isFree()) { value.bg -= factor * h[row] * Y[row];  ++row; }
        if (value.cc.isFree()) { value.cc -= factor * h[row] * Y[row];  ++row; }
        
        for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
          if (value.parameters[pmt].QE .isFree()) { value.parameters[pmt].QE  -= factor * h[row] * Y[row];  ++row; }
          if (value.parameters[pmt].TTS.isFree()) { value.parameters[pmt].TTS -= factor * h[row] * Y[row];  ++row; }
          if (value.parameters[pmt].t0 .isFree()) { value.parameters[pmt].t0  -= factor * h[row] * Y[row];  ++row; }
        }
      }

      return { precessor / estimator->getRho(1.0), ndf };
    }

void JCALIBRATE::JFit::evaluate ( const data_type &  data )

inlineprivate

Evaluation of fit.

Parameters

data

data

Definition at line 1330 of file JFitK40.hh.

    {
      using namespace std;
      using namespace JPP;
      
      typedef JModel::real_type real_type;


      successor = 0.0;

      V.reset();
      Y.reset();


      // model parameter indices

      const struct M_t {
        M_t(const JModel& model)
        {
          R  = model.getIndex(&JK40Parameters_t::R);
          p1 = model.getIndex(&JK40Parameters_t::p1);
          p2 = model.getIndex(&JK40Parameters_t::p2);
          p3 = model.getIndex(&JK40Parameters_t::p3);
          p4 = model.getIndex(&JK40Parameters_t::p4);
          bg = model.getIndex(&JK40Parameters_t::bg);
          cc = model.getIndex(&JK40Parameters_t::cc);
        }

        int R;
        int p1;
        int p2;
        int p3;
        int p4;
        int bg;
        int cc;

      } M(value);


      // PMT parameter indices

      struct I_t {
        I_t(const JModel& model, const int pmt) :
          QE (INVALID_INDEX),
          TTS(INVALID_INDEX),
          t0 (INVALID_INDEX)
        {
          const int index = model.getIndex(pmt);

          int N = 0;

          if (model.parameters[pmt].QE .isFree()) { QE  = index + N; ++N; }
          if (model.parameters[pmt].TTS.isFree()) { TTS = index + N; ++N; }
          if (model.parameters[pmt].t0 .isFree()) { t0  = index + N; ++N; }
        }

        int QE;
        int TTS;
        int t0;
      };


      typedef vector< pair<int, double> >  buffer_type;
 
      buffer_type  buffer;

      for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {

        const pair_type& pair = ix->first;

        if (value.parameters[pair.first ].status &&
            value.parameters[pair.second].status) {

          const real_type& real = value.getReal(pair);

          const JGauss gauss(real.t0, real.sigma, real.signal);

          const double            R1  = value.getValue   (real.ct);
          const JK40Parameters_t& R1p = value.getGradient(real.ct);

          const std::pair<I_t, I_t> PMT(I_t(value, pair.first),
                                        I_t(value, pair.second));

          for (const rate_type& iy : ix->second) {

            const double  R2  = gauss.getValue   (iy.dt_ns);
            const JGauss& R2p = gauss.getGradient(iy.dt_ns);        

            const double  R   = value.bg()  +  R1 * (value.cc() + R2);
            const double  u   =  (iy.value - R)       / iy.error;
            const double  w   = -estimator->getPsi(u) / iy.error;

            successor += estimator->getRho(u);

            buffer.clear();

            if (M.R  != INVALID_INDEX) { buffer.push_back({M.R,  w * (value.cc() + R2) * R1p.R () * value.R .getDerivative()}); }
            if (M.p1 != INVALID_INDEX) { buffer.push_back({M.p1, w * (value.cc() + R2) * R1p.p1() * value.p1.getDerivative()}); }
            if (M.p2 != INVALID_INDEX) { buffer.push_back({M.p2, w * (value.cc() + R2) * R1p.p2() * value.p2.getDerivative()}); }
            if (M.p3 != INVALID_INDEX) { buffer.push_back({M.p3, w * (value.cc() + R2) * R1p.p3() * value.p3.getDerivative()}); }
            if (M.p4 != INVALID_INDEX) { buffer.push_back({M.p4, w * (value.cc() + R2) * R1p.p4() * value.p4.getDerivative()}); }
            if (M.bg != INVALID_INDEX) { buffer.push_back({M.bg, w *  value.bg.getDerivative()}); }
            if (M.cc != INVALID_INDEX) { buffer.push_back({M.cc, w *  R1 * R1p.cc() * value.cc.getDerivative()}); }

            if (PMT.first .QE  != INVALID_INDEX) { buffer.push_back({PMT.first .QE , w * R1 * R2p.signal * value.parameters[pair.second].QE () * value.parameters[pair.first ].QE .getDerivative()}); }
            if (PMT.second.QE  != INVALID_INDEX) { buffer.push_back({PMT.second.QE , w * R1 * R2p.signal * value.parameters[pair.first ].QE () * value.parameters[pair.second].QE .getDerivative()}); }
            if (PMT.first .TTS != INVALID_INDEX) { buffer.push_back({PMT.first .TTS, w * R1 * R2p.sigma  * value.parameters[pair.first ].TTS() * value.parameters[pair.first ].TTS.getDerivative() / real.sigma}); }
            if (PMT.second.TTS != INVALID_INDEX) { buffer.push_back({PMT.second.TTS, w * R1 * R2p.sigma  * value.parameters[pair.second].TTS() * value.parameters[pair.second].TTS.getDerivative() / real.sigma}); }
            if (PMT.first .t0  != INVALID_INDEX) { buffer.push_back({PMT.first .t0,  w * R1 * R2p.mean   * value.parameters[pair.first ].t0 .getDerivative() * +1.0}); }
            if (PMT.second.t0  != INVALID_INDEX) { buffer.push_back({PMT.second.t0,  w * R1 * R2p.mean   * value.parameters[pair.second].t0 .getDerivative() * -1.0}); }

            for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {

              Y[row->first] += row->second;

              V[row->first][row->first] += row->second * row->second;

              for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
                V[row->first][col->first] += row->second * col->second;
                V[col->first][row->first]  = V[row->first][col->first];
              }
            }
          }
        }
      }
    }

Member Data Documentation

constexpr int JCALIBRATE::JFit::MAXIMUM_ITERATIONS = 10000

static

maximal number of iterations.

Definition at line 1308 of file JFitK40.hh.

constexpr double JCALIBRATE::JFit::EPSILON = 1.0e-4

static

maximal distance to minimum.

Definition at line 1309 of file JFitK40.hh.

constexpr double JCALIBRATE::JFit::LAMBDA_MIN = 0.01

static

minimal value control parameter

Definition at line 1310 of file JFitK40.hh.

constexpr double JCALIBRATE::JFit::LAMBDA_MAX = 100.0

static

maximal value control parameter

Definition at line 1311 of file JFitK40.hh.

constexpr double JCALIBRATE::JFit::LAMBDA_UP = 10.0

static

multiplication factor control parameter

Definition at line 1312 of file JFitK40.hh.

constexpr double JCALIBRATE::JFit::LAMBDA_DOWN = 10.0

static

multiplication factor control parameter

Definition at line 1313 of file JFitK40.hh.

constexpr double JCALIBRATE::JFit::PIVOT = std::numeric_limits<double>::epsilon()

static

minimal value diagonal element of matrix

Definition at line 1314 of file JFitK40.hh.

int JCALIBRATE::JFit::debug

Definition at line 1316 of file JFitK40.hh.

estimator_type JCALIBRATE::JFit::estimator

M-Estimator function.

Definition at line 1317 of file JFitK40.hh.

double JCALIBRATE::JFit::lambda

Definition at line 1319 of file JFitK40.hh.

JModel JCALIBRATE::JFit::value

Definition at line 1320 of file JFitK40.hh.

int JCALIBRATE::JFit::numberOfIterations

Definition at line 1321 of file JFitK40.hh.

JMATH::JMatrixNS JCALIBRATE::JFit::V

Definition at line 1322 of file JFitK40.hh.

JMATH::JVectorND JCALIBRATE::JFit::Y

private

Definition at line 1457 of file JFitK40.hh.

double JCALIBRATE::JFit::successor

private

Definition at line 1458 of file JFitK40.hh.

JModel JCALIBRATE::JFit::previous

private

Definition at line 1459 of file JFitK40.hh.

std::vector<double> JCALIBRATE::JFit::h

private

Definition at line 1460 of file JFitK40.hh.

---

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

• JFitK40.hh