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

#include <JKatoomba.hh>

Inheritance diagram for JACOUSTICS::JKatoomba< JGandalf >:

Public Types
typedef double	result_type

Public Member Functions
	JKatoomba (const JDetector &detector, const JSoundVelocity &velocity, const int option)
	Constructor. More...

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

Public Attributes
double	lambda

JModel	value

int	numberOfIterations

JMATH::JMatrixNS	V

Static Public Attributes
static int	debug = 0
	debug level More...

static int	MAXIMUM_ITERATIONS = 1000
	maximal number of iterations More...

static double	EPSILON = 1.0e-3
	maximal distance to minimum More...

static double	LAMBDA_MIN = 0.01
	minimal value control parameter More...

static double	LAMBDA_MAX = 100.0
	maximal value control parameter More...

static double	LAMBDA_UP = 9.0
	multiplication factor control parameter More...

static double	LAMBDA_DOWN = 11.0
	multiplication factor control parameter More...

static double	PIVOT = 1.0e-3
	minimal value diagonal element of matrix More...

Private Member Functions
template<class T >
void	evaluate (T __begin, T __end)
	Evaluation of fit. More...

Private Attributes
std::vector< double >	Y

result_type	successor

JModel	previous

std::vector< double >	h

Detailed Description

template<>
struct JACOUSTICS::JKatoomba< JGandalf >

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

Definition at line 552 of file JKatoomba.hh.

Member Typedef Documentation

typedef double JACOUSTICS::JKatoomba< JGandalf >::result_type

Definition at line 555 of file JKatoomba.hh.

Constructor & Destructor Documentation

JACOUSTICS::JKatoomba< JGandalf >::JKatoomba	(	const JDetector &	detector,
		const JSoundVelocity &	velocity,
		const int	option
	)

inline

Constructor.

Parameters

detector	detector
velocity	sound velocity
option	option

Definition at line 565 of file JKatoomba.hh.

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

Member Function Documentation

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 580 of file JKatoomba.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);
 
       lambda    = LAMBDA_MIN;
 
       result_type precessor = numeric_limits<double>::max();
       
       for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
 
         DEBUG("step:     " << numberOfIterations << endl);
 
         evaluate(__begin, __end);
 
         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(__begin, __end);
         }
 
 
         // 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;
         }
 
 
         try {
           V.invert();
         }
         catch (const JException& error) {
           ERROR("JKatoomb<JGandalf>: " << error.what() << endl);
           return precessor;
         }
 
 
         for (int i = 0; i != N; ++i) {
 
           DEBUG("u[" << noshowpos << setw(3) << i << "] = " << showpos << FIXED(20,5) << value[i]);
 
           double y = 0.0;
 
           for (int j = 0; j != N; ++j) {
             y += V(i,j) * Y[j] * h[i] * h[j];
           }
 
           value[i] -= y;
 
           DEBUG(' ' << FIXED(20,10) << y << noshowpos << endl);
         }
       }
 
       return precessor;
     }

template<class T >

void JACOUSTICS::JKatoomba< JGandalf >::evaluate	(	T	__begin,
		T	__end
	)

inlineprivate

Evaluation of fit.

Parameters

__begin	begin of data
__end	end of data

Definition at line 711 of file JKatoomba.hh.

     {
       using namespace std;
       using namespace JPP;
 
       successor = 0.0;
 
       V.reset();
 
       for (std::vector<double>::iterator i = Y.begin(); i != Y.end(); ++i) {
         *i = 0.0;
       }
 
       for (T hit = __begin; hit != __end; ++hit) {
 
         const JGEOMETRY::JString& string     = detector[hit->getString()];
         const JMODEL   ::JString& parameters = value.string[hit->getString()];
         const JPosition3D         position   = string.getPosition(parameters, hit->getFloor());
 
         const double D     = hit->getDistance(position);
         const double Vi    = velocity.getInverseVelocity(D, hit->getZ(), position.getZ());
         const double toa_s = value.emitter[hit->getEKey()].t1  +  D * Vi;
           
         const double u     = (toa_s - hit->getValue()) / hit->sigma;
         const double W     = sqrt(hit->getWeight());
 
         successor += (W*W) * estimator->getRho(u);
         
         H_t H(1.0, string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);
 
         H    *= W * estimator->getPsi(u) / hit->sigma;
 
         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)  += H.vs  * H.t1;
             V(i.vs,  i.tx)  += H.vs  * H.tx;
             V(i.vs,  i.ty)  += H.vs  * H.ty;
             V(i.vs,  i.tx2) += H.vs  * H.tx2;
             V(i.vs,  i.ty2) += H.vs  * H.ty2;
             
             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)  += H.tx2 * H.t1;
             V(i.tx2, i.tx)  += H.tx2 * H.tx;
             V(i.tx2, i.ty)  += H.tx2 * H.ty;
 
             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)  += H.ty2 * H.t1;
             V(i.ty2, i.tx)  += H.ty2 * H.tx;
             V(i.ty2, i.ty)  += H.ty2 * H.ty;
 
             V(i.tx2, i.tx2) += H.tx2 * H.tx2;  V(i.tx2, i.ty2) += H.tx2 * H.ty2;
             V(i.ty2, i.tx2) += H.ty2 * H.tx2;  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)  += H.tx  * H.t1;   V(i.ty,  i.t1)  += H.ty  * H.t1;
 
             V(i.tx,  i.tx)  += H.tx  * H.tx;   V(i.tx,  i.ty)  += H.tx  * H.ty;
             V(i.ty,  i.tx)  += H.ty  * H.tx;   V(i.ty,  i.ty)  += H.ty  * H.ty;
 
             Y[i.tx]  += W * H.tx;
             Y[i.ty]  += W * H.ty;
             break;
           
           default:
             break;
           }
         }
       }
     }