Template specialisation of fit function of acoustic model based on JGandalf minimiser. More...

#include <JKatoomba_t.hh>

Inheritance diagram for JACOUSTICS::JKatoomba< JGandalf >:

Public Types
typedef double	result_type

typedef std::map< JLocation, std::map< JEmitter, std::vector< JHit > > >	data_type
	Type definition internal data structure.

Public Member Functions
	JKatoomba (const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
	Constructor The option corresponds to the use of fit parameters in the model of the detector geometry.

template<class T >
result_type	operator() (T __begin, T __end)
	Fit.

Public Attributes
double	lambda

JModel	value

int	numberOfIterations

JMATH::JMatrixNS	V

Static Public Attributes
static int	debug
	debug level

static int	MAXIMUM_ITERATIONS
	maximal number of iterations

static double	EPSILON
	maximal distance to minimum

static double	LAMBDA_MIN
	minimal value control parameter

static double	LAMBDA_MAX
	maximal value control parameter

static double	LAMBDA_UP
	multiplication factor control parameter

static double	LAMBDA_DOWN
	multiplication factor control parameter

static double	PIVOT
	minimal value diagonal element of matrix

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

Private Attributes
JMATH::JVectorND	Y

result_type	successor

JModel	previous

std::vector< double >	h

Detailed Description

Template specialisation of fit function of acoustic model based on JGandalf minimiser.

Definition at line 678 of file JKatoomba_t.hh.

Member Typedef Documentation

◆ result_type

double JACOUSTICS::JKatoomba< JGandalf >::result_type

Definition at line 681 of file JKatoomba_t.hh.

◆ data_type

std::map<JLocation, std::map<JEmitter, std::vector<JHit> > > JACOUSTICS::JKatoomba< JGandalf >::data_type

Type definition internal data structure.

Definition at line 686 of file JKatoomba_t.hh.

Constructor & Destructor Documentation

◆ JKatoomba()

JACOUSTICS::JKatoomba< JGandalf >::JKatoomba	(	const JGeometry &	geometry,
		const JSoundVelocity &	velocity,
		const int	option )

inline

Constructor The option corresponds to the use of fit parameters in the model of the detector geometry.

A negative implies that all strings in the detector use common fit parameters.

Parameters

geometry	detector geometry
velocity	sound velocity
option	option

Definition at line 698 of file JKatoomba_t.hh.

                                            :
      JKatoomba<>(geometry, velocity, option)
    {};

Member Function Documentation

◆ operator()()

template<class T >

result_type JACOUSTICS::JKatoomba< JGandalf >::operator()	(	T	__begin,
		T	__end )

inline

Fit.

Parameters

__begin	begin of hits
__end	end of hits

Returns: chi2 and gradient

Definition at line 713 of file JKatoomba_t.hh.

    {
      using namespace std;
      using namespace JPP;
 
      value.setOption(this->option);
 
      const int N = value.getN();
 
      V.resize(N);
      Y.resize(N);
      h.resize(N);
 
      data_type data;
 
      for (T hit = __begin; hit != __end; ++hit) {
        data[hit->getLocation()][hit->getEmitter()].push_back(*hit);
      }
 
      lambda    = LAMBDA_MIN;
 
      result_type 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,5) << successor << endl);
 
        if (successor < precessor) {
 
          if (numberOfIterations != 0) {
 
            if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
              return successor;
            }
 
            if (lambda > LAMBDA_MIN) {
              lambda /= LAMBDA_DOWN;
            }
          }
 
          precessor = successor;
          previous  = value;
 
        } else {
 
          value   = previous;
          lambda *= LAMBDA_UP;
 
          if (lambda > LAMBDA_MAX) {
            return precessor;                                   // no improvement found
          }
 
          evaluate(data);
        }
 
        // force definite positiveness
 
        for (int i = 0; i != N; ++i) {
 
          if (V(i,i) < PIVOT) {
            V(i,i) = PIVOT;
          }
 
          h[i] = 1.0 / sqrt(V(i,i));
        }
 
        // normalisation
 
        for (int i = 0; i != N; ++i) {
          for (int j = 0; j != i; ++j) {
            V(j,i) *= h[i] * h[j];
            V(i,j)  = V(j,i);
          }
        }
 
        for (int i = 0; i != N; ++i) {
          V(i,i) = 1.0 + lambda;
        }
 
 
        // solve A x = b
 
        for (int 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
 
        for (int row = 0; row != N; ++row) {
 
          DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(20,5) << value[row]);
 
          value[row] -= h[row] * Y[row];
 
          DEBUG(" -> " << FIXED(20,10) << value[row] << noshowpos << endl);
        }
      }
 
      return precessor;
    }

◆ evaluate()

void JACOUSTICS::JKatoomba< JGandalf >::evaluate ( const data_type & data )

inlineprivate

Evaluation of fit.

Parameters

data data

Definition at line 849 of file JKatoomba_t.hh.

    {
      using namespace std;
      using namespace JPP;
 
      successor = 0.0;
 
      V.reset();
      Y.reset();
 
      for (data_type::const_iterator p = data.begin(); p != data.end(); ++p) {
 
        const JGEOMETRY::JString& string     = geometry    [p->first.getString()];
        const JMODEL   ::JString& parameters = value.string[p->first.getString()];
        const JPosition3D         position   = string.getPosition(parameters, p->first.getFloor());
 
        for (data_type::mapped_type::const_iterator emitter = p->second.begin(); emitter != p->second.end(); ++emitter) {
 
          const double D  = emitter->first.getDistance(position);
          const double Vi = velocity.getInverseVelocity(D, emitter->first.getZ(), position.getZ());
 
          const H_t    H0(1.0, string.getGradient(parameters, emitter->first.getPosition(), p->first.getFloor()) * Vi);
 
          for (data_type::mapped_type::mapped_type::const_iterator hit = emitter->second.begin(); hit != emitter->second.end(); ++hit) {
        
            const double toa_s = value.emission[hit->getEKey()].t1  +  D * Vi;
          
            const double u  = (toa_s - hit->getValue()) / hit->getSigma();
            const double W  = sqrt(hit->getWeight());
 
            successor += (W*W) * estimator->getRho(u);
        
            const H_t    H  = H0 * (W * estimator->getPsi(u) / hit->getSigma());
 
            I_t i;
 
            i.t1  =  value.getIndex(hit->getEKey(),   &H_t::t1);
            i.tx  =  value.getIndex(hit->getString(), &H_t::tx);
            i.ty  =  value.getIndex(hit->getString(), &H_t::ty);
            i.tx2 =  value.getIndex(hit->getString(), &H_t::tx2);
            i.ty2 =  value.getIndex(hit->getString(), &H_t::ty2);
            i.vs  =  value.getIndex(hit->getString(), &H_t::vs);
 
            V(i.t1, i.t1) += H.t1 * H.t1;
 
            Y[i.t1] += W * H.t1;
 
            if (hit->getFloor() != 0) {
 
              switch (this->option) {
 
              case JMODEL::FIT_EMITTERS_AND_STRINGS_2nd_ORDER_AND_STRETCHING_t:
                V(i.t1,  i.vs)  += H.t1  * H.vs;  V(i.tx,  i.vs)  += H.tx  * H.vs;  V(i.ty,  i.vs)  += H.ty  * H.vs;  V(i.tx2, i.vs)  += H.tx2 * H.vs;  V(i.ty2, i.vs)  += H.ty2 * H.vs;
 
                V(i.vs,  i.t1)   = V(i.t1,  i.vs);
                V(i.vs,  i.tx)   = V(i.tx,  i.vs);
                V(i.vs,  i.ty)   = V(i.ty,  i.vs);
                V(i.vs,  i.tx2)  = V(i.tx2, i.vs);
                V(i.vs,  i.ty2)  = V(i.ty2, i.vs);
 
                V(i.vs,  i.vs)  += H.vs  * H.vs;
 
                Y[i.vs]  += W * H.vs;
            
              case JMODEL::FIT_EMITTERS_AND_STRINGS_2nd_ORDER_t:
                V(i.t1,  i.tx2) += H.t1  * H.tx2;  V(i.tx,  i.tx2) += H.tx  * H.tx2;  V(i.ty,  i.tx2) += H.ty  * H.tx2;
 
                V(i.tx2, i.t1)   = V(i.t1,  i.tx2);
                V(i.tx2, i.tx)   = V(i.tx,  i.tx2);
                V(i.tx2, i.ty)   = V(i.ty,  i.tx2);
            
                V(i.t1,  i.ty2) += H.t1  * H.ty2;  V(i.tx,  i.ty2) += H.tx  * H.ty2;  V(i.ty,  i.ty2) += H.ty  * H.ty2;
 
                V(i.ty2, i.t1)   = V(i.t1,  i.ty2);
                V(i.ty2, i.tx)   = V(i.tx,  i.ty2);
                V(i.ty2, i.ty)   = V(i.ty,  i.ty2);
 
                V(i.tx2, i.tx2) += H.tx2 * H.tx2;    V(i.tx2, i.ty2) += H.tx2 * H.ty2;
                V(i.ty2, i.tx2)  = V(i.tx2, i.ty2);  V(i.ty2, i.ty2) += H.ty2 * H.ty2;
 
                Y[i.tx2] += W * H.tx2;
                Y[i.ty2] += W * H.ty2;
 
              case JMODEL::FIT_EMITTERS_AND_STRINGS_1st_ORDER_t:
                V(i.t1,  i.tx)  += H.t1  * H.tx;   V(i.t1,  i.ty)  += H.t1  * H.ty;
                V(i.tx,  i.t1)   = V(i.t1, i.tx);  V(i.ty,  i.t1)   = V(i.t1, i.ty);
 
                V(i.tx,  i.tx)  += H.tx  * H.tx;   V(i.tx,  i.ty)  += H.tx  * H.ty;
                V(i.ty,  i.tx)   = V(i.tx, i.ty);  V(i.ty,  i.ty)  += H.ty  * H.ty;
 
                Y[i.tx]  += W * H.tx;
                Y[i.ty]  += W * H.ty;
                break;
          
              default:
                break;
              }
            }
          }
        }
      }
    }