/*
 * Deriv.hh
 *
 *  Created on: Dec 3, 2012
 *      Author: f.casimir
 */

#ifndef Deriv_HH_
#define Deriv_HH_

#include "Parameter.hh"
#include "ParamData.hh"
#include "DataTypeMath.hh"
#include "Operation.hh"

namespace AMDA {
namespace Parameters {

/**
 * @class Deriv
 * @brief It is responsible to compute the deriv data.
 * @details This class implement the interface Operation.
 * @code
 * 		3 point Lagrangian interpolation (points {x0,y0; x1,y1; x2,y2}):
 * 		df/dx = y0*(2x-x1-x2)/(x01*x02) + y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1)/(x20*x21),
 * 		where
 * 		x01= x0-x1; x02=x0-x2, x12=x1-x2...
 * @endcode
 */
template<class TParamData>
class Deriv : public Operation {
public:
  /**
   * @brief  Element type of paramData
   */
  typedef typename TParamData::ElementType ElementType;

  /**
   * @brief Constructor.
   * @details Create the ParamData type of the input ParamData.
   */
  Deriv(Process& pProcess, TParamData& paramInput)
	: Operation(pProcess),
	  _paramInput(paramInput),
	  _paramOutput(new TParamData()),
	  _firstPointComputed(false) {
	  _paramDataOutput=_paramOutput;
  }

  /**
   * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
   */
  void  write(ParamDataIndexInfo &pParamDataIndexInfo) {

	//Store data in a working list
	for (unsigned int index = pParamDataIndexInfo._startIndex;
		index< pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess; index++) {
		_data.push_back(std::make_pair(_paramInput.getTime(index), _paramInput.getDataList()[index]));
	}

	ElementType res;
	typename std::list<std::pair<double, ElementType> >::iterator itData0, itData1, itData2;

	if (!_firstPointComputed && (_data.size() >= 3)) {
		//First point calculation
		itData0 = _data.begin();
		itData1 = std::next(itData0, 1);
		itData2 = std::next(itData0, 2);
		res = computeDeriv(*itData0, *itData1, *itData2);
		_paramOutput->getDataList().push_back(res);
		_paramOutput->pushTime((*itData0).first);
		_firstPointComputed = true;
	}

	while (_data.size() >= 3) {
		//Calculation for other points than the first one or the last one
		itData0 = _data.begin();
		itData1 = std::next(itData0, 1);
		itData2 = std::next(itData0, 2);
		res = computeDeriv(*itData0, *itData1, *itData2);
		_paramOutput->getDataList().push_back(res);
		_paramOutput->pushTime((*itData1).first);
		_lastPopData = *itData0;
		_data.pop_front();
	}

	if (pParamDataIndexInfo._timeIntToProcessChanged || pParamDataIndexInfo._noMoreTimeInt) {
		//Finish process for this interval
		if (!_firstPointComputed) {
			//Only one or two points => cannot compute deriv => set NaN values
			while (!_data.empty()) {
				res << NotANumber();
				_paramOutput->getDataList().push_back(res);
				_paramOutput->pushTime((*_data.begin()).first);
				_data.pop_front();
			}

		}
		else {
			//Last point calculation
			itData1 = _data.begin();
			itData2 = std::next(itData1, 1);
			res = computeDeriv(_lastPopData, *itData1, *itData2);
			_paramOutput->getDataList().push_back(res);
			_paramOutput->pushTime((*itData2).first);
			_data.clear();
		}
	}
  }

  ElementType computeDeriv(std::pair<double, ElementType> data0, std::pair<double, ElementType> data1, std::pair<double, ElementType> data2) {
	ElementType y0 = data0.second;
	ElementType y1 = data1.second;
	ElementType y2 = data2.second;
	double x0 = data0.first;
	double x1 = data1.first;
	double x2 = data2.first;

	double x = x1;
	double x02 = x0 -x2;
	double x01 = x0 -x1;
	double x10 = x1 -x0;
	double x12 = x1 -x2;
	double x20 =x2 - x0 ;
	double x21 = x2 -x1;
	double _2x = 2*x;

	return y0*(_2x-x1-x2)/(x01*x02) + y1*(_2x-x0-x2)/(x10*x12)+y2*(_2x-x0-x1)/(x20*x21);
  }
  
  /**
    * @overload Operation::reset(double pStartTime, double pTimeInt)
    * @brief reset static data to process another TimeInterval
    */
  virtual  void  reset() {
        Operation::reset();
	_data.clear();
        _firstPointComputed = false;
  }

private:
	/**<
	 * @brief It is the channel of data derived
	 */
  TParamData &_paramInput;
	/**<
	 * @brief It is the channel of the data derived
	 */
  TParamData *_paramOutput;

  bool _firstPointComputed;
  std::pair<double, ElementType> _lastPopData;
  std::list<std::pair<double, ElementType> > _data;
};

} /* namespace Parameters */
} /* namespace AMDA */
#endif /* Deriv_HH_ */