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 "JTools/JHashMapToolkit.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_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
    };

    inline JOption_t getOption();
    inline void setOption(const int);

    /**
     * Auxiliary class to hide access function to fit option.
     */
    struct JModel_t {

      friend JOption_t JMODEL::getOption();
      friend void JMODEL::setOption(const int);

    private:
      /**
       * Get reference to fit option.
       *
       * \return                    option
       */
      static JOption_t& get_option()
      {
        static JOption_t option = fit_emitters_and_strings_1st_order_t;

        return option;
      }
    };
    

    /**
     * Get fit option.
     *
     * \return                    option
     */
    inline JOption_t getOption()
    {
      return JModel_t::get_option();
    }


    /**
     * Set fit option.
     *
     * \param  option             option
     */
    inline void setOption(const int option)
    {
      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:
        JModel_t::get_option() = static_cast<JOption_t>(option);
        break;

      default:
        THROW(JValueOutOfRange, "Invalid option " << option);
      }
    }


    /**
     * 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.
       *
       * \return                          number of parameters
       */
      static inline size_t getN()
      {
        switch (JMODEL::getOption()) {

        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 " << JMODEL::getOption());
        }
      }


      /**
       * 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;
    };


    /**
     * Emitter parameters.
     */
    struct JEmitter :
      public JMath  <JEmitter>,
      public JEquals<JEmitter>
    {
      /**
       * Default constructor.
       */
      JEmitter() :
        t1(0.0)
      {}


      /**
       * Constructor.
       *
       * \param  t1                       time-of-emission [s]
       */
      JEmitter(const double t1) :
        t1(t1)
      {}


      /**
       * Get number of fit parameters.
       *
       * \return                          number of parameters
       */
      static inline size_t getN()
      {
        return 1;
      }


      /**
       * Negate emitter.
       *
       * \return                          this emitter
       */
      JEmitter& negate()
      {
        t1 = -t1;

        return *this;
      }


      /**
       * Add emitter.
       *
       * \param  emitter                  emitter
       * \return                          this emitter
       */
      JEmitter& add(const JEmitter& emitter)
      {
        t1 += emitter.t1;

        return *this;
      }


      /**
       * Subtract emitter.
       *
       * \param  emitter                  emitter
       * \return                          this emitter
       */
      JEmitter& sub(const JEmitter& emitter)
      {
        t1 -= emitter.t1;

        return *this;
      }


      /**
       * Scale emitter.
       *
       * \param  factor                   multiplication factor
       * \return                          this emitter
       */
      JEmitter& mul(const double factor)
      {
        t1 *= factor;

        return *this;
      }


      /**
       * Scale emitter.
       *
       * \param  factor                   division factor
       * \return                          this emitter
       */
      JEmitter& div(const double factor)
      {
        t1 /= factor;

        return *this;
      }


      /**
       * Check equality.
       *
       * \param  emitter                  emitter
       * \param  precision                precision
       * \return                          true if emitters are equal; else false
       */
      bool equals(const JEmitter& emitter,
                  const double    precision = std::numeric_limits<double>::min()) const
      {
        return (fabs(t1 - emitter.t1) <= precision);
      }


      /**
       * Write emitter parameters to output stream.
       *
       * \param  out                      output stream
       * \param  emitter                  emitter
       * \return                          output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const JEmitter& emitter)
      {
        using namespace std;
        using namespace JPP;

        return out << FIXED(20,6) << emitter.t1;
      }


      double t1;
    };
  }


  /**
   * Model for fit to acoustics data.
   *
   * The model consists of string parameters and emitter parameters.\n
   * These parameters relate to the string identifer and emitter key, respectively.
   */
  struct JModel :
    public JMath  <JModel>,
    public JEquals<JModel>,
    public JTOOLS::JHashMapToolkit
  {
    using JTOOLS::JHashMapToolkit::reset;
    using JTOOLS::JHashMapToolkit::set;
    using JTOOLS::JHashMapToolkit::equals;

    typedef JMODEL::JString   JString;
    typedef JMODEL::JEmitter  JEmitter;


    /**
     * Type definition of fit parameter.
     */
    typedef size_t  parameter_type;


    /**
     * Default constructor.
     */
    JModel()
    {}


    /**
     * 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 emitter key
     * </pre>
     *
     * \param  __begin                  begin of hits
     * \param  __end                    end   of hits
     */
    template<class T>
    JModel(T __begin, T __end)
    {
      for (T hit = __begin; hit != __end; ++hit) {

        if (!this->string.has(hit->getString())) {
          this->string[hit->getString()] = JString();
        }

        if (!this->emitter.has(hit->getEKey())) {
          this->emitter[hit->getEKey()]  = JEmitter();
        }
      }
    }


    /**
     * Reset parameters.
     *
     * \param  zero                     zero
     * \return                          this model
     */
    JModel& operator=(const JMATH::JZero& zero)
    {
      this->reset();

      return *this;
    }


    /**
     * Clear parameters.
     */
    void clear()
    {
      string .clear();
      emitter.clear();
    }


    /**
     * Reset parameters.
     */
    void reset()
    {
      reset(string);
      reset(emitter);
    }


    /**
     * Negate model.
     *
     * \return                          this model
     */
    JModel& negate()
    {
      evaluate(this->string,  &JString ::negate);
      evaluate(this->emitter, &JEmitter::negate);

      return *this;
    }


    /**
     * Add model.
     *
     * \param  model                    model
     * \return                          this model
     */
    JModel& add(const JModel& model)
    {
      evaluate(this->string,  model.string,  &JString ::add);
      evaluate(this->emitter, model.emitter, &JEmitter::add);

      return *this;
    }
    

    /**
     * Subtract model.
     *
     * \param  model                    model
     * \return                          this model
     */
    JModel& sub(const JModel& model)
    {
      evaluate(this->string,  model.string,  &JString ::sub);
      evaluate(this->emitter, model.emitter, &JEmitter::sub);

      return *this;
    }


    /**
     * Scale model.
     *
     * \param  factor                   multiplication factor
     * \return                          this model
     */
    JModel& mul(const double factor)
    {
      evaluate(this->string,  &JString ::mul, factor);
      evaluate(this->emitter, &JEmitter::mul, factor);

      return *this;
    }


    /**
     * Scale model.
     *
     * \param  factor                   division factor
     * \return                          this model
     */
    JModel& div(const double factor)
    {
      evaluate(this->string,  &JString ::div, factor);
      evaluate(this->emitter, &JEmitter::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 (equals(this->string,  model.string,  precision) &&
              equals(this->emitter, model.emitter, 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 (JHashMap<int, JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
        out << "string:  " << setw(4) << i->first << ' ' << i->second << endl;
      }

      for (JHashMap<JEKey, JEmitter>::const_iterator i = model.emitter.begin(); i != model.emitter.end(); ++i) {
        out << "emitter: " << setw(3) << i->first << ' ' << i->second << endl;
      }

      return out;
    }


    /**
     * Get number of fit parameters.
     *
     * \return                          number of parameters
     */
    size_t getN() const
    {
      return emitter.getN() + string.getN();
    }


    /**
     * 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
    {
      if (string.has(id))
        return emitter.getN()  +  string.getIndex(id) * JString::getN()  +  getIndex(p);
      else
        THROW(JValueOutOfRange, "Invalid identifier " << id);
    }


    /**
     * Get index of fit parameter for given emitter.
     *
     * \param  id                       emitter key
     * \param  p                        pointer to data member
     * \return                          parameter
     */
    size_t getIndex(const JEKey& id, double JEmitter::*p) const
    {
      if (emitter.has(id))
        return emitter.getIndex(id) * JEmitter::getN()  +  getIndex(p);
      else
        THROW(JValueOutOfRange, "Invalid identifier " << id);
    }


    /**
     * 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 < emitter.getN()) { return getParameter(emitter, i); }

      i -= emitter.getN();

      if (i < string .getN()) { return getParameter(string,  i); }

      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 < emitter.getN()) { return getParameter(emitter, i); }

      i -= emitter.getN();

      if (i < string .getN()) { return getParameter(string,  i); }

      THROW(JIndexOutOfRange, "Invalid index " << index << " >= " << getN());
    }


    /**
     * Map emitter key to model parameters of emitter.
     */
    struct emitter_type :
      public JHashMap<JEKey, JEmitter>
    {
    private:
      /**
       * Auxiliary class for multiple associative map operators.
       */
      struct JHashMapHelper {
        /**
         * Constructor.
         *
         * \param  map                      map
         * \param  id                       emitter identifier
         */
        JHashMapHelper(JHashMap<JEKey, JEmitter>& map, int id) :
          map(map),
          id (id)
        {}


        /**
         * Get value corresponding to event counter (i.e.\ second part of JEKey).
         *
         * \param  counter                  event counter
         * \return                          emitter
         */
        JEmitter& operator[](const int counter)
        {
          return map[JEKey(id, counter)];
        }

      private:
        JHashMap<JEKey, JEmitter>& map;
        const int id;
      };

    public:

      using JHashMap<JEKey, JEmitter>::operator[];

      /**
       * Get number of fit parameters.
       *
       * \return                          number of parameters
       */
      size_t getN() const
      {
        return this->size() * JEmitter::getN();
      }

      
      /**
       * Get helper corresponding to emitter identifier (i.e.\ first part of JEKey).
       *
       * \param  id                       emitter identifier
       * \return                          helper
       */
      JHashMapHelper operator[](int id) 
      {
        return JHashMapHelper(*this, id);
      }

    } emitter;


    /**
     * Map string identifier to model parameters of string.
     */
    struct string_type :
      public JHashMap<int, JString>
    {
      /**
       * Get number of fit parameters.
       *
       * \return                          number of parameters
       */
      size_t getN() const
      {
        return this->size() * JString::getN();
      }
    } string;

  private:
    /**
     * Get index of fit parameter in given data structure.
     *
     * \param  p                        pointer to data member
     * \return                          index
     */
    template<class T>
    static inline size_t getIndex(double T::*p) 
    {
      T* __p__ = NULL;

      return ((double*) &(__p__->*p) - (double*) __p__);
    }


    /**
     * Get read access to fit parameter value at given index in buffer.
     *
     * \param  buffer                   buffer
     * \param  index                    index
     * \return                          value
     */
    template<class JKey_t, class JValue_t, class JEvaluator_t>
    static double getParameter(const JHashMap<JKey_t, JValue_t, JEvaluator_t>& buffer, const size_t index)
    {
      const size_t pos       = index / JValue_t::getN();         // position of element in buffer
      const size_t offset    = index % JValue_t::getN();         // offset of parameter in element

      const JValue_t& value  = buffer.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  buffer                   buffer
     * \param  index                    index
     * \return                          value
     */
    template<class JKey_t, class JValue_t, class JEvaluator_t>
    static double& getParameter(JHashMap<JKey_t, JValue_t, JEvaluator_t>& buffer, const size_t index)
    {
      const size_t pos       = index / JValue_t::getN();         // position of element in buffer
      const size_t offset    = index % JValue_t::getN();         // offset of parameter in element

      JValue_t&    value     = buffer.data()[pos].second;        // value of element at given position

      return ((double*) &value)[offset];                         // parameter at given offset
    }
  };
}

#endif