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. More...

Public Member Functions
	JKatoomba (const JGeometry &geometry, 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 = 10.0
	multiplication factor control parameter More...

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

static 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
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 735 of file JKatoomba_t.hh.

Member Typedef Documentation

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

Definition at line 738 of file JKatoomba_t.hh.

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

Type definition internal data structure.

Definition at line 743 of file JKatoomba_t.hh.

Constructor & Destructor Documentation

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

inline

Constructor.

Parameters

geometry	detector geometry
velocity	sound velocity
option	option

Definition at line 753 of file JKatoomba_t.hh.

                                             :
       JKatoomba<>(geometry, 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 768 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;
     }

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

inlineprivate

Evaluation of fit.

Parameters

data data

Definition at line 904 of file JKatoomba_t.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 (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;
               }
             }
           }
         }
       }
     }