1 #ifndef __JACOUSTICS__JKATOOMBA__
2 #define __JACOUSTICS__JKATOOMBA__
29 namespace JACOUSTICS {}
30 namespace JPP {
using namespace JACOUSTICS; }
32 namespace JACOUSTICS {
77 template<
class JPDF_t>
126 static bool option =
true;
217 template<
template<
class T>
class JMinimiser_t>
253 template<
class JPDF_t>
256 const double toa_s = this->getToA(model, hit);
257 const double u = (toa_s - hit.getValue()) / hit.sigma;
259 return estimator->getRho(u);
313 using namespace JGEOMETRY;
315 value =
JModel(__begin, __end);
320 const size_t N = (
getOption() ? value.getN() : value.emitter.getN());
326 for (
T hit = __begin; hit != __end; ++hit) {
330 const double Vi = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.
getZ());
332 const double h1 =
string.getHeight(hit->getFloor());
333 const JPosition3D p1 =
string.getPosition() - hit->getPosition();
335 const double y = hit->getValue() - Vi*ds;
336 const double W = sqrt(hit->getWeight());
339 H.tx = W * Vi * p1.
getX() * h1 / ds;
340 H.ty = W * Vi * p1.
getY() * h1 / ds;
342 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
343 i.tx = value.getIndex(hit->getString(), &H_t::tx);
344 i.ty = value.getIndex(hit->getString(), &H_t::ty);
346 V(i.t1, i.t1) +=
H.t1 *
H.t1;
348 Y[i.t1] += W *
H.t1 * y;
352 V(i.t1, i.tx) +=
H.t1 *
H.tx; V(i.t1, i.ty) +=
H.t1 *
H.ty;
353 V(i.tx, i.t1) +=
H.tx *
H.t1; V(i.ty, i.t1) +=
H.ty *
H.t1;
355 V(i.tx, i.tx) +=
H.tx *
H.tx; V(i.tx, i.ty) +=
H.tx *
H.ty;
356 V(i.ty, i.tx) +=
H.ty *
H.tx; V(i.ty, i.ty) +=
H.ty *
H.ty;
358 Y[i.tx] += W *
H.tx * y;
359 Y[i.ty] += W *
H.ty * y;
367 for (
size_t row = 0; row !=
N; ++row) {
368 for (
size_t col = 0; col !=
N; ++col) {
369 value[row] += V(row,col) * Y[col];
409 template<
class JPDF_t>
412 const double toa_s = this->getToA(model, hit);
413 const double u = (toa_s - hit.getValue()) / hit.sigma;
415 return estimator->getRho(u);
437 this->step.push_back(model);
441 this->step.push_back(model);
450 this->step.push_back(model);
497 const int N = (
getOption() ? value.getN() : value.emitter.getN());
505 result_type precessor = numeric_limits<double>::max();
507 for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
509 DEBUG(
"step: " << numberOfIterations << endl);
511 evaluate(__begin, __end);
513 DEBUG(
"lambda: " <<
FIXED(12,5) << lambda << endl);
514 DEBUG(
"chi2: " <<
FIXED(12,5) << successor << endl);
516 if (successor < precessor) {
518 if (numberOfIterations != 0) {
520 if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
524 if (lambda > LAMBDA_MIN) {
525 lambda /= LAMBDA_DOWN;
529 precessor = successor;
537 if (lambda > LAMBDA_MAX) {
541 evaluate(__begin, __end);
547 for (
int i = 0; i !=
N; ++i) {
549 if (V(i,i) < PIVOT) {
553 h[i] = 1.0 / sqrt(V(i,i));
559 for (
int i = 0; i !=
N; ++i) {
560 for (
int j = 0;
j != i; ++
j) {
561 V(
j,i) *= h[i] * h[
j];
566 for (
int i = 0; i !=
N; ++i) {
567 V(i,i) = 1.0 + lambda;
575 ERROR(
"JKatoomb<JGandalf>: " << error.
what() << endl);
580 for (
int i = 0; i !=
N; ++i) {
582 DEBUG(
"u[" << noshowpos << i <<
"] = " << showpos <<
FIXED(15,5) << value[i]);
586 for (
int j = 0;
j !=
N; ++
j) {
587 y += V(i,
j) *
Y[
j] * h[i] * h[
j];
592 DEBUG(
' ' <<
FIXED(15,5) << y << noshowpos << endl);
634 for (
T hit = __begin; hit != __end; ++hit) {
640 const double D = hit->getDistance(position);
641 const double Vi = velocity.getInverseVelocity(D, hit->getZ(), position.
getZ());
642 const double toa_s = value.emitter[hit->getEKey()].t1 + D * Vi;
644 const double u = (toa_s - hit->getValue()) / hit->sigma;
645 const double W = sqrt(hit->getWeight());
647 successor += (W*W) * estimator->getRho(u);
649 H_t H(1.0,
string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);
651 H *= W * estimator->getPsi(u) / hit->sigma;
655 i.
t1 = value.getIndex(hit->getEKey(), &H_t::t1);
656 i.
tx = value.getIndex(hit->getString(), &H_t::tx);
657 i.
ty = value.getIndex(hit->getString(), &H_t::ty);
Auxiliary base class for fit function of acoustic model.
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity)
Constructor.
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity)
Constructor.
JFit_t
Enumeration for fit algorithms.
General purpose data regression method.
JString & mul(const double factor)
Scale string.
Wrapper class around STL string class.
Auxiliary base class for aritmetic operations of derived class types.
int getFloor() const
Get floor number.
static double LAMBDA_UP
multiplication factor control parameter
Interface for maximum likelihood estimator (M-estimator).
result_type operator()(const JFunction_t &fit, T __begin, T __end)
Get chi2.
JParserTemplateElement< bool > getOption(bool &object, const std::string &name, const std::string &help="")
Auxiliary method for creation of template parser element object.
I_t(const int t1, const int tx, const int ty)
Constructor.
Indices of H-equation in global model.
static const double H
Planck constant [eV s].
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
JEmitter & mul(const double factor)
Scale emitter.
then for HISTOGRAM in h0 h1
double operator()(const JModel &model, const JHit< JPDF_t > &hit) const
Fit function.
I_t()
Default constructor.
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Auxiliary data structure for floating point format specification.
JEKey getEKey() const
Get emitter hash key of this hit.
Template definition of linear fit.
double getDistance(const JVector3D &pos) const
Get distance to point.
Model for fit to acoustics data.
H_t(const JMODEL::JEmitter &emitter, const JMODEL::JString &string)
Constructor.
static double EPSILON
maximal distance to minimum
The template JSharedPointer class can be used to share a pointer to an object.
result_type operator()(T __begin, T __end)
Fit.
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
JKatoomba_t(const JDetector &detector, const JSoundVelocity &velocity)
Constructor.
static double LAMBDA_MIN
minimal value control parameter
void evaluate(T __begin, T __end)
Evaluation of fit.
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity)
Constructor.
virtual double getInverseVelocity(const double D_m, const double z1, const double z2) const override
Get inverse velocity of sound.
static int debug
debug level
do set_variable OUTPUT_DIRECTORY $WORKDIR T
double operator()(T __begin, T __end)
Fit.
result_type operator()(const JModel &model, const JHit< JPDF_t > &hit) const
Fit function.
double getY() const
Get y position.
const JPosition3D & getPosition() const
Get position.
static bool & get_option()
Get reference to fit option.
double getLengthSquared() const
Get length squared.
General purpose messaging.
Implementation for velocity of sound.
static double LAMBDA_MAX
maximal value control parameter
static int MAXIMUM_ITERATIONS
maximal number of iterations
double operator()(const JFunction_t &fit, T __begin, T __end)
Multi-dimensional fit.
void invert()
Invert matrix according LDU decomposition.
JACOUSTICS::JModel::string_type string
JACOUSTICS::JModel::emitter_type emitter
int getString() const
Get string number.
virtual const char * what() const override
Get error message.
const JModel & operator()(T __begin, T __end) const
Get start values of string parameters.
static double LAMBDA_DOWN
multiplication factor control parameter
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++...
result_type operator()(const JModel &model, T __begin, T __end)
Fit.
H_t()
Default constructor.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
static double PIVOT
minimal value diagonal element of matrix
double getToA(const JModel &model, const JHit< JPDF_t > &hit) const
Get estimated time-of-arrival for given hit.
double getX() const
Get x position.
H_t & mul(const double factor)
Scale H-equation.
static void setOption(const bool option)
Set fit option.
Data structure for position in three dimensions.
const JDetector & detector
Model for fit to acoutsics data.
static bool getOption()
Get fit option.
then usage $script[input file[working directory[option]]] nWhere option can be N
do echo Generating $dir eval D
double getZ() const
Get z position.
Maximum likelihood estimator (M-estimators).
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity)
Constructor.
Template definition of fit function of acoustic model.