#ifndef __JMATH__JMODEL__
#define __JMATH__JMODEL__

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

#include "JLang/JEquals.hh"
#include "JLang/JManip.hh"
#include "JMath/JMath.hh"
#include "JMath/JZero.hh"

/**
 * \author mdejong
 */

namespace JMATH {}
namespace JPP { using namespace JMATH; }

namespace JMATH {

  using JLANG::JEquals;

  /**
   * Fit model.
   */
  struct JModel_t :
    public std::vector<double>,
    JMath<JModel_t>
  {
    using std::vector<double>::operator[];

    typedef size_t      parameter_type;

    /**
     * Reset.
     *
     * \param  zero          zero
     * \return               this model
     */
    JModel_t& operator=(const JMATH::JZero& zero)
    {
      for (size_t i = 0; i != this->size(); ++i) {
	(*this)[i] = 0.0;
      }

      return *this;
    }


    /**
     * Get value at given index.
     *
     * \param  index         index
     * \return               value
     */
    double& operator[](const size_t index) 
    {
      if (index >= this->size()) {
	this->resize(index + 1, 0.0);
      }

      return std::vector<double>::operator[](index);
    }


    /**
     * Equality.
     *
     * \param  model         model
     * \param  eps           numerical precision
     * \return               true if model's identical; else false
     */
    bool equals(const JModel_t& model,
                const double    eps = std::numeric_limits<double>::min()) const
    {
      for (size_t i = 0; i != this->size(); ++i) {
	if (fabs((*this)[i] - model[i]) >= eps) {
	  return false;
	}
      }

      return true;
    }


    /**
     * Negate model.
     *
     * \return               this model
     */
    JModel_t& negate()
    {
      for (size_t i = 0; i != this->size(); ++i) {
	(*this)[i] = -(*this)[i];
      }

      return *this;
    }


    /**
     * Add model.
     *
     * \param  model        model
     * \return               this model
     */
    JModel_t& add(const JModel_t& model)
    {
      for (size_t i = 0; i != model.size(); ++i) {
	(*this)[i] += model[i];
      }

      return *this;
    }


    /**
     * Subtract model.
     *
     * \param  model        model
     * \return               this model
     */
    JModel_t& sub(const JModel_t& model)
    {
      for (size_t i = 0; i != model.size(); ++i) {
	(*this)[i] -= model[i];
      }

      return *this;
    }


    /**
     * Scale model.
     *
     * \param  factor        multiplication factor
     * \return               this model
     */
    JModel_t& mul(const double factor)
    {
      for (size_t i = 0; i != this->size(); ++i) {
	(*this)[i] *= factor;
      }

      return *this;
    }


    /**
     * Scale model.
     *
     * \param  factor        division factor
     * \return               this model
     */
    JModel_t& div(const double factor)
    {
      for (size_t i = 0; i != this->size(); ++i) {
	(*this)[i] /= factor;
      }

      return *this;
    }


    /**
     * Write model to input stream.
     *
     * \param  in            input stream
     * \param  model         model
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JModel_t& model)
    {
      model.clear();

      for (double value; in >> value; ) {
	model.push_back(value);
      }

      return in;
    }


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

      for (JModel_t::const_iterator i = model.begin(); i != model.end(); ++i) {
	out << FIXED(7,3) << *i  << ' ';
      }

      return out;
    }
  };
}

#endif