#ifndef __JACOUSTICS__JMODEL__ #define __JACOUSTICS__JMODEL__ #include #include #include #include #include #include "JMath/JMath.hh" #include "JMath/JZero.hh" #include "JLang/JEquals.hh" #include "JLang/JException.hh" #include "JLang/JManip.hh" #include "JTools/JHashMap.hh" #include "JAcoustics/JEKey.hh" /** * \file * * Model for fit to acoutsics data. * \author mdejong */ namespace JACOUSTICS {} namespace JPP { using namespace JACOUSTICS; } namespace JACOUSTICS { using JMATH::JMath; using JLANG::JEquals; using JLANG::JValueOutOfRange; using JLANG::JIndexOutOfRange; using JTOOLS::JHashMap; /** * Auxiliary namespace to encapsulate different model parameters.\n */ namespace JMODEL { /** * Fit options. */ enum JOption_t { FIT_UNDEFINED_t = -1, //!< fit undefined FIT_EMITTERS_ONLY_t = 0, //!< fit only times of emission of emitters FIT_EMITTERS_AND_STRINGS_1st_ORDER_t = 1, //!< fit times of emission of emitters and tilt angles of strings FIT_EMITTERS_AND_STRINGS_2nd_ORDER_t = 2, //!< fit times of emission of emitters and tilt angles of strings with second order correction FIT_EMITTERS_AND_STRINGS_2nd_ORDER_AND_STRETCHING_t = 3 //!< fit times of emission of emitters and tilt angles of strings with second order correction and stretching of string }; /** * String parameters. */ struct JString : public JMath , public JEquals { /** * Default constructor. */ JString() : tx (0.0), ty (0.0), tx2(0.0), ty2(0.0), vs (0.0) {} /** * Constructor. * * \param tx slope dx/dz * \param ty slope dy/dz * \param tx2 2nd order correction of slope dx/dz * \param ty2 2nd order correction of slope dy/dz * \param vs stretching factor */ JString(const double tx, const double ty, const double tx2 = 0.0, const double ty2 = 0.0, const double vs = 0.0) : tx (tx), ty (ty), tx2(tx2), ty2(ty2), vs (vs) {} /** * Get number of fit parameters. * * \param option option * \return number of parameters */ static inline size_t getN(const JMODEL::JOption_t option) { switch (option) { case FIT_UNDEFINED_t: case FIT_EMITTERS_ONLY_t: return 0; case FIT_EMITTERS_AND_STRINGS_1st_ORDER_t: return 2; case FIT_EMITTERS_AND_STRINGS_2nd_ORDER_t: return 4; case FIT_EMITTERS_AND_STRINGS_2nd_ORDER_AND_STRETCHING_t: return 5; default: THROW(JValueOutOfRange, "Invalid option " << option); } } /** * Negate string. * * \return this string */ JString& negate() { tx = -tx; ty = -ty; tx2 = -tx2; ty2 = -ty2; vs = -vs; return *this; } /** * Add string. * * \param string string * \return this string */ JString& add(const JString& string) { tx += string.tx; ty += string.ty; tx2 += string.tx2; ty2 += string.ty2; vs += string.vs; return *this; } /** * Subtract string. * * \param string string * \return this string */ JString& sub(const JString& string) { tx -= string.tx; ty -= string.ty; tx2 -= string.tx2; ty2 -= string.ty2; vs -= string.vs; return *this; } /** * Scale string. * * \param factor multiplication factor * \return this string */ JString& mul(const double factor) { tx *= factor; ty *= factor; tx2 *= factor; ty2 *= factor; vs *= factor; return *this; } /** * Scale string. * * \param factor division factor * \return this string */ JString& div(const double factor) { tx /= factor; ty /= factor; tx2 /= factor; ty2 /= factor; vs /= factor; return *this; } /** * Check equality. * * \param string string * \param precision precision * \return true if strings are equal; else false */ bool equals(const JString& string, const double precision = std::numeric_limits::min()) const { return (fabs(tx - string.tx) <= precision && fabs(ty - string.ty) <= precision && fabs(tx2 - string.tx2) <= precision && fabs(ty2 - string.ty2) <= precision && fabs(vs - string.vs) <= precision); } /** * Get length squared. * * \return square of length */ double getLengthSquared() const { return tx*tx + ty*ty; } /** * Get length. * * \return length */ double getLength() const { return sqrt(getLengthSquared()); } /** * Get angle. * * \return angle [rad] */ double getAngle() const { return atan2(ty, tx); } /** * Get dot product. * * \param string string * \return dot product */ double getDot(const JString& string) const { return (tx * string.tx + ty * string.ty + tx2 * string.tx2 + ty2 * string.ty2 + vs * string.vs); } /** * Read string parameters from input stream. * * \param in input stream * \param string string * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JString& string) { return in >> string.tx >> string.ty >> string.tx2 >> string.ty2 >> string.vs; } /** * Write string parameters to output stream. * * \param out output stream * \param string string * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const JString& string) { using namespace std; using namespace JPP; return out << FIXED(10,7) << string.tx << ' ' << FIXED(10,7) << string.ty << ' ' << SCIENTIFIC(12,3) << string.tx2 << ' ' << SCIENTIFIC(12,3) << string.ty2 << ' ' << FIXED(8,5) << string.vs; } double tx; double ty; double tx2; double ty2; double vs; }; /** * Emission parameters. */ struct JEmission : public JMath , public JEquals { /** * Default constructor. */ JEmission() : t1(0.0) {} /** * Constructor. * * \param t1 time-of-emission [s] */ JEmission(const double t1) : t1(t1) {} /** * Get number of fit parameters. * * \param option option * \return number of parameters */ static inline size_t getN(const JMODEL::JOption_t option) { switch (option) { case FIT_UNDEFINED_t: return 0; default: return 1; } } /** * Negate emission. * * \return this emission */ JEmission& negate() { t1 = -t1; return *this; } /** * Add emission. * * \param emission emission * \return this emission */ JEmission& add(const JEmission& emission) { t1 += emission.t1; return *this; } /** * Subtract emission. * * \param emission emission * \return this emission */ JEmission& sub(const JEmission& emission) { t1 -= emission.t1; return *this; } /** * Scale emission. * * \param factor multiplication factor * \return this emission */ JEmission& mul(const double factor) { t1 *= factor; return *this; } /** * Scale emission. * * \param factor division factor * \return this emission */ JEmission& div(const double factor) { t1 /= factor; return *this; } /** * Check equality. * * \param emission emission * \param precision precision * \return true if emissions are equal; else false */ bool equals(const JEmission& emission, const double precision = std::numeric_limits::min()) const { return (fabs(t1 - emission.t1) <= precision); } /** * Write emission parameters to output stream. * * \param out output stream * \param emission emission * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const JEmission& emission) { using namespace std; using namespace JPP; return out << FIXED(20,6) << emission.t1; } double t1; }; } /** * Model for fit to acoustics data. * * The model consists of string parameters and emission parameters.\n * These parameters relate to the string identifer and emission key, respectively. */ struct JModel : public JMath , public JEquals { typedef JMODEL::JString JString; typedef JMODEL::JEmission JEmission; /** * Type definition of fit parameter. */ typedef size_t parameter_type; /** * Default constructor. */ JModel() : option(JMODEL::FIT_UNDEFINED_t) {} /** * Constructor. * * This constructor can be used to set up a default model (i.e. all values at zero) for the given set of hits.\n * The data type corresponding to the hits should provide for the following policy methods. *
     *    int      getString();        // get string identifier
     *    JEKey    getEKey();          // get emission key
     * 
* * \param __begin begin of hits * \param __end end of hits */ template JModel(T __begin, T __end) : option(JMODEL::FIT_UNDEFINED_t) { for (T hit = __begin; hit != __end; ++hit) { if (!this->string.has(hit->getString())) { this->string[hit->getString()] = JString(); } if (!this->emission.has(hit->getEKey())) { this->emission[hit->getEKey()] = JEmission(); } } } /** * Reset parameters. * * \param zero zero * \return this model */ JModel& operator=(const JMATH::JZero& zero) { this->reset(); return *this; } /** * Get fit option. * * \return option */ JMODEL::JOption_t getOption() const { return option; } /** * Set fit option. * * \param option option */ void setOption(const int option) { using namespace JMODEL; switch (option) { case FIT_EMITTERS_ONLY_t: case FIT_EMITTERS_AND_STRINGS_1st_ORDER_t: case FIT_EMITTERS_AND_STRINGS_2nd_ORDER_t: case FIT_EMITTERS_AND_STRINGS_2nd_ORDER_AND_STRETCHING_t: this->option = static_cast(option); break; default: THROW(JValueOutOfRange, "Invalid option " << option); } } /** * Clear parameters. */ void clear() { string .clear(); emission.clear(); } /** * Reset parameters. */ void reset() { string .reset(); emission.reset(); } /** * Negate model. * * \return this model */ JModel& negate() { this->string .evaluate(&JString ::negate); this->emission.evaluate(&JEmission::negate); return *this; } /** * Add model. * * \param model model * \return this model */ JModel& add(const JModel& model) { this->string .evaluate(model.string, &JString ::add); this->emission.evaluate(model.emission, &JEmission::add); return *this; } /** * Subtract model. * * \param model model * \return this model */ JModel& sub(const JModel& model) { this->string .evaluate(model.string, &JString ::sub); this->emission.evaluate(model.emission, &JEmission::sub); return *this; } /** * Scale model. * * \param factor multiplication factor * \return this model */ JModel& mul(const double factor) { this->string .evaluate(&JString ::mul, factor); this->emission.evaluate(&JEmission::mul, factor); return *this; } /** * Scale model. * * \param factor division factor * \return this model */ JModel& div(const double factor) { this->string .evaluate(&JString ::div, factor); this->emission.evaluate(&JEmission::div, factor); return *this; } /** * Check equality. * * \param model model * \param precision precision * \return true if models are equal; else false */ bool equals(const JModel& model, const double precision = std::numeric_limits::min()) const { return (this->string .equals(model.string, precision) && this->emission.equals(model.emission, precision)); } /** * Write model parameters to output stream. * * \param out output stream * \param model model * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const JModel& model) { using namespace std; using namespace JPP; for (hash_map::const_iterator i = model.string.begin(); i != model.string.end(); ++i) { out << "string: " << setw(4) << i->first << ' ' << i->second << endl; } for (hash_map::const_iterator i = model.emission.begin(); i != model.emission.end(); ++i) { out << "emission: " << setw(3) << i->first << ' ' << i->second << endl; } return out; } /** * Get number of fit parameters. * * \return number of parameters */ size_t getN() const { return emission.getN(this->option) + string.getN(this->option); } /** * Get index of fit parameter for given string. * * \param id string identifier * \param p pointer to data member * \return parameter */ size_t getIndex(int id, double JString::*p) const { return emission.getN(this->option) + string.getIndex(id, p, this->option); } /** * Get index of fit parameter for given emission. * * \param id emission key * \param p pointer to data member * \return parameter */ size_t getIndex(const JEKey& id, double JEmission::*p) const { return emission.getIndex(id, p, this->option); } /** * Read access to fit parameter value by index. * * \param index index * \return value */ inline double operator[](const size_t index) const { size_t i = index; if (i < emission.getN(this->option)) { return emission.getParameter(i, this->option); } i -= emission.getN(this->option); if (i < string .getN(this->option)) { return string .getParameter(i, this->option); } THROW(JIndexOutOfRange, "Invalid index " << index << " >= " << getN()); } /** * Read/write access to fit parameter value by index. * * \param index index * \return value */ inline double& operator[](const size_t index) { size_t i = index; if (i < emission.getN(this->option)) { return emission.getParameter(i, this->option); } i -= emission.getN(this->option); if (i < string .getN(this->option)) { return string .getParameter(i, this->option); } THROW(JIndexOutOfRange, "Invalid index " << index << " >= " << getN()); } /** * Auxiliary data structure with extended functionality of hash-map. */ template struct hash_map : public JHashMap { /** * Get number of fit parameters. * * \param option option * \return number of parameters */ size_t getN(const JMODEL::JOption_t option) const { return this->size() * value_type::getN(option); } /** * Get index of parameter. * * \param key key * \param p pointer to data member * \param option option * \return index */ size_t getIndex(const key_type key, double value_type::*p, const JMODEL::JOption_t option) const { value_type* __p__ = NULL; if (this->has(key)) return static_cast&>(*this).getIndex(key) * value_type::getN(option) + ((double*) &(__p__->*p) - (double*) __p__); else THROW(JValueOutOfRange, "Invalid key " << key); } /** * Get read access to fit parameter value at given index in buffer. * * \param index index * \param option option * \return value */ double getParameter(const size_t index, const JMODEL::JOption_t option) const { const size_t pos = index / value_type::getN(option); // position of element in buffer const size_t offset = index % value_type::getN(option); // offset of parameter in element const value_type& value = this->data()[pos].second; // value of element at given position return ((const double*) &value)[offset]; // parameter at given offset } /** * Get read/write access to fit parameter value at given index in buffer. * * \param index index * \param option option * \return value */ double& getParameter(const size_t index, const JMODEL::JOption_t option) { const size_t pos = index / value_type::getN(option); // position of element in buffer const size_t offset = index % value_type::getN(option); // offset of parameter in element value_type& value = this->data()[pos].second; // value of element at given position return ((double*) &value)[offset]; // parameter at given offset } /** * Evaluate arithmetic operation. * * \param f1 operation */ void evaluate(value_type& (value_type::*f1)()) { for (typename JHashMap::iterator i = this->begin(); i != this->end(); ++i) { (i->second.*f1)(); } } /** * Evaluate arithmetic operation. * * \param buffer buffer * \param f1 operation */ void evaluate(const hash_map& buffer, value_type& (value_type::*f1)(const value_type&)) { for (typename JHashMap::const_iterator i = buffer.begin(); i != buffer.end(); ++i) { ((*this)[i->first].*f1)(i->second); } } /** * Evaluate arithmetic operation. * * \param f1 operation * \param factor factor */ void evaluate(value_type& (value_type::*f1)(const double), const double factor) { for (typename JHashMap::iterator i = this->begin(); i != this->end(); ++i) { (i->second.*f1)(factor); } } /** * Check equality of hash map. * * \param buffer buffer * \return true if hash map is equal; else false */ bool equals(const hash_map& buffer) const { for (typename JHashMap::const_iterator p = this->begin(), q = buffer.begin(); ; ++p, ++q) { if (p != this->end() && q != buffer.end()) { if (p->first != q->first || p->second != q->second) { return false; } } else if (p == this->end() && q == buffer.end()) { return true; } else { return false; } } } /** * Check equality of has map. * * \param buffer buffer * \param precision precision * \return true if hash map is equal; else false */ bool equals(const hash_map& buffer, const double precision) const { for (typename JHashMap::const_iterator p = this->begin(), q = buffer.begin(); ; ++p, ++q) { if (p != this->end() && q != buffer.end()) { if (p->first != q->first || !p->second.equals(q->second, precision)) { return false; } } else if (p == this->end() && q == buffer.end()) { return true; } else { return false; } } } }; /** * Map emission key to model parameters of emission. */ struct emission_type : public hash_map { private: /** * Auxiliary class for multiple associative map operators. */ struct helper { /** * Constructor. * * \param map map * \param id emission identifier */ helper(hash_map& map, int id) : map(map), id (id) {} /** * Get value corresponding to event counter (i.e.\ second part of JEKey). * * \param counter event counter * \return emission */ JEmission& operator[](const int counter) { return map[JEKey(id, counter)]; } private: hash_map& map; const int id; }; public: using hash_map::operator[]; /** * Get helper corresponding to emission identifier (i.e.\ first part of JEKey). * * \param id emission identifier * \return helper */ helper operator[](int id) { return helper(*this, id); } } emission; /** * Map string identifier to model parameters of string. */ struct string_type : public hash_map {} string; private: JMODEL::JOption_t option; }; } #endif