JAcoustics/JModel.hh

Go to the documentation of this file.

#ifndef __JACOUSTICS__JMODEL__
#define __JACOUSTICS__JMODEL__

#include <istream>
#include <ostream>
#include <iomanip>
#include <limits>
#include <cmath>

#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;
  using JTOOLS::JHashEvaluator;


  /**
   * 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  <JString>,
      public JEquals<JString>
    {
      /**
       * 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<double>::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  <JEmission>,
      public JEquals<JEmission>
    {
      /**
       * 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<double>::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  <JModel>,
    public JEquals<JModel>
  {
    typedef JMODEL::JString   JString;
    typedef JMODEL::JEmission JEmission;


    /**
     * 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.
     * <pre>
     *    int      getString();        // get string identifier
     *    JEKey    getEKey();          // get emission key
     * </pre>
     *
     * \param  __begin                  begin of hits
     * \param  __end                    end   of hits
     */
    template<class T>
    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<JOption_t>(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<double>::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<int, JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
        out << "string:   " << setw(4) << i->first << ' ' << i->second << endl;
      }

      for (hash_map<JEKey, JEmission>::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 for common fit parameters.
     */
    struct hash_evaluator {

      /**
       * Option for common fit parameters.
       */
      static bool singularity;
      
      inline int operator()(const char key) const                   { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const unsigned char key) const          { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const short key) const                  { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const unsigned short key) const         { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const int key) const                    { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const unsigned int key) const           { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const long int key) const               { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const unsigned long int key) const      { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const long long int key) const          { return (singularity ? 0 : (int) key); }  //!< Get hash value.
      inline int operator()(const unsigned long long int key) const { return (singularity ? 0 : (int) key); }  //!< Get hash value.
    };
    

    /**
     * Auxiliary data structure with extended functionality of hash-map.
     */
    template<class key_type, class value_type, class evaluator_type = JHashEvaluator>
    struct hash_map :
      public JHashMap<key_type, value_type, evaluator_type>
    {
      typedef JHashMap<key_type, value_type, evaluator_type>  hashmap_type;

      /**
       * 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<const hashmap_type&>(*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 hashmap_type::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 hashmap_type::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 hashmap_type::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 hashmap_type::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 hashmap_type::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<JEKey, JEmission>
    {
    private:
      /**
       * Auxiliary class for multiple associative map operators.
       */
      struct helper {
        /**
         * Constructor.
         *
         * \param  map                      map
         * \param  id                       emission identifier
         */
        helper(hash_map<JEKey, JEmission>& 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<JEKey, JEmission>& map;
        const int id;
      };

    public:

      using hash_map<JEKey, JEmission>::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<int, JString, hash_evaluator>
    {} string;

  private:
    JMODEL::JOption_t option;
  };
}

#endif