//------------------------------------------------------------------------------
//	Copyright (C) 2003, 2004, 2009 Humboldt-Universitaet zu Berlin
//
//	This library is free software; you can redistribute it and/or
//	modify it under the terms of the GNU Lesser General Public
//	License as published by the Free Software Foundation; either
//	version 2.1 of the License, or (at your option) any later version.
//
//	This library is distributed in the hope that it will be useful,
//	but WITHOUT ANY WARRANTY; without even the implied warranty of
//	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//	Lesser General Public License for more details.
//
//	You should have received a copy of the GNU Lesser General Public
//	License along with this library; if not, write to the Free Software
//	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//------------------------------------------------------------------------------

/**
 * @file base/Continuous.h
 * @author Ralf Gerstenberger
 * @date created at 2003/06/15
 * @brief Declaration of odemx::base::Continuous, its observer, and odemx::base::ContuTrace
 * @sa Continuous.cpp
 * @since 1.0
 */

#ifndef ODEMX_BASE_CONTINUOUS_INCLUDED
#define ODEMX_BASE_CONTINUOUS_INCLUDED

#include <odemx/setup.h>

// usage of this class requires SimTime type double,
// which can be switched at compile time by defining
// ODEMX_USE_CONTINUOUS in file odemx/setup.h
#ifdef ODEMX_USE_CONTINUOUS

#include <odemx/base/Process.h>
#include <odemx/base/TypeDefs.h>

#include <string>
#include <vector>
#include <list>

namespace odemx {
namespace base {

// forward declaration
class ContinuousObserver;

/** \class Continuous

	\ingroup base

	\author Ralf Gerstenberger

	\brief Time continuous process

	\note %Continuous supports Observation.
	\note %Continuous supports Trace.

	\sa Process

	%Continuous is a base class for all user defined time continuous
	processes. It provides functionality of Process and introduces
	the possibility to define time continuous processes. A user has
	to implement the derivatives() function to define the time
	continuous state changes. The function integrate() computes
	the time continuous state changes.

	\par Example:
	\include continuousProcesses.cpp

	\since 1.0
*/
class Continuous
:	public Process // is DataProducer, etc.
,	public data::Observable< ContinuousObserver >
{
public:

	/**
		\brief Construction for user-defined Simulation
		\param s
			pointer to the Simulation object
		\param l
			label of this object
		\param d
			number of state variables used in derivatives()
		\param o
			initial observer
	*/
	Continuous( Simulation& sim, const data::Label& label, int dimension,
			ContinuousObserver* o = 0 );

	/// Destruction
	virtual ~Continuous();

	/**
		\brief Interrupt

		If a %Continuous process is interrupted inside integrate()
		the computation is not synchronised to the time of the
		interrupt. There might be an error depending on the step length.
		This can be prevented if the process that caused the
		interrupt is registered as a peer to the continuous process.

		\sa Process::interrupt
	*/
	virtual void interrupt();

	/**
		\brief Get number of state variables

		%getDimension() returns the number of state variables specified
		in the constructor. No more than the specified state variables
		may me used in derivatives(). The state variables are stored in
		the STL vector state.
	*/
	unsigned int getDimension() const;

	/**
		\brief Set minimum and maximum step length

		\param min
			minimum step length
		\param max
			maximum step length

		The integrate() function computes state changes step by step.
		The step length used is variable depending on numerical	errors,
		state events and registered peer processes. With %setStepLength()
		the user can choose boundaries for the step length according to
		the needed accuracy.

		\sa integrate(), addPeer()

	*/
	void setStepLength( double min, double max );

	/**
		\brief Get current step length

		The integrate() function computes state changes step by step.
		The step length used is variable depending on numerical	errors,
		state events and registered peer processes. %getStepLength()
		returns the current step length.

		 \sa integrate()
	*/
	double getStepLength() const;

	/**
		\brief Set errorlimit

		\param type
			error type: 1 for relative error; 0 for absolute error
		\param limit
			error limit

		The integrate() function computes state changes step by step.
		The step length used does depend on the numerical error during
		computation. With %setErrorlimit() the user can set the highest
		acceptable error. If the actual error is higher the step length
		is reduced. \n
		The limit can be set absolute or relative. In the second case the
		limit is interpreted relative to the actual state values.

		\sa integrate()

	*/
	void setErrorlimit( int type, double limit );

	/**
		\brief Add a partner process

		\param p
			pointer to Process object

		The integrate() function computes state changes step by step. The
		next state value is computed in advance to the time consumption.
		If another process is accessing the state variables during the
		time consumption it receives wrong values. The step length however
		is variable and can be changed to synchronise computation with
		partner Processes. The function  %addPeer() is used to register
		partner Processes for synchronisation.

		\note
		A peer must have a higher priority than its	continuous partner
		to interact properly. To support this a continuous process
		reduces its priority by one inside integrate(). A peer has to
		ensure that its priority is at least as high as the priority of
		the continuous process.

		\sa integrate(), delPeer()
	*/
	void addPeer( Process* newPeer );

	/**
		\brief Remove a peer

		\param p
			pointer to Process object

		Removes a registered partner process.

		\sa addPeer()
	*/
	void delPeer( Process* newPeer );

	/**
		\brief State variables

		The \c state variables used during time continuous state changes
		are stored in this STL vector. A process which is accessing
		this variables directly should be registered with addPeer().

		\sa integrate(), derivatives(), addPeer()
	*/
	std::vector< double > state;

	/**
		\brief Rate variables

		The STL vector \c rate is used for the computation	of time continuous
		state changes. The state variables are not changed directly. In
		derivatives() the user sets the rate by witch each state variable
		is changed.

		\sa integrate(), derivatives()
	*/
	std::vector< double > rate;

	/**
		\brief Error vector

		The step by step computation of time continuous state changes
		introduces an inherent error. This error is computed for every state
		variable independently. \c error_vector is used to store the error
		values. The virtual function errorNorm() is called to compute a
		single error value from this vector which is stored in \c error_value.

		\sa integrate(), errorNorm(), error_value
	*/
	std::vector< double > error_vector;

	/**
		\brief Error

		The step by step computation of time continuous state changes
		introduces an inherent error. This error is computed for every state
		variable independently. The function errorNorm() is used to compute
		a single error value from the vector of error values which is stored
		in \c error_value.

		\sa integrate(), errorNorm()
	*/
	double error_value;

protected:
	/**
		\brief User defined process behaviour

		This function is implemented by the user to	define the process
		behaviour.

		\sa Process::main()
	*/
	virtual int main() = 0;

	/**
		\brief User defined continuous behaviour

		\param t
			time

		The user has to implement the function %derivatives() to define
		the time continuous state changes. This is done indirectly by
		setting the rate by witch each state variable changes. The function
		is called by integrate() several times for each step.

		\dontinclude continuousExample.cpp
		\par Example Oscillator
		The %derivatives() function describes a simple sin/cos oscillator.
		The two rate variables depend only on the two state variables:
		\skipline virtual void derivatives (double t) {
		\until }

		\par Example FreeFall
		This %derivatives() function uses the provided parameter \p t
		to include the current time in the computation:
		\skipline virtual void derivatives (double t) {
		\until }

		\note
		If the time is required for the computation of the rates the
		prameter \p t has to be used instead of getCurrentTime().

		\sa integrate(), rate, state
	*/
    virtual void derivatives( double t ) = 0;

	/**
		\brief Computation of time continuous state changes

		\param timeEvent
			time limit for integration (0->infinite)
		\param stateEvent
			state event that stops integration by returning true (member call-back-function)
		\return 0 if timeEvent hit; 1 if stateEvent hit;
				2 if interrupted

		The function %integrate() computes time continuous state changes
		using the Runge-Kutta algorithm. The computation is stopped when
		the time \p timeEvent is reached, the state event occurred or
		the interrupt() function was called. \n
		The computation is done step by step with a varying step length. The
		actual step length depends on computation errors, registered peers
		and the stop time (\p timeEvent). The user can set boundaries to
		the step length with the function setStepLength(). The default
		values used by integrate are 0.01 and 0.1 for minimum and maximum.
		\n
		The step by step computation of time continuous processes introduces
		an inherent error. integrate() observes the computational errors for
		each state variable. The function errorNorm() is called to compute
		a single value from all error values. The user can use the function
		setErrorlimit() to specify the largest acceptable error. If the single
		value computed by errorNorm() exceeds the acceptable error %integrate()
		tries to reduce the step length. This fails if the resulting step
		length would be smaller than the minimum step length set with setStepLength().
		The default error limit is an absolute error of 0.1.
		\n
		If during the computation a state event occurred, integrate() starts
		a binary search to find the exact time of the event. The minimum
		step length set by setStepLength() defines the accuracy of this
		search.
		\n
		%integrate() synchronises computation with partner processes. A process
		is registered as a partner with the function addPeer().

		\sa derivatives(), setStepLength(), setErrorlimit(), errorNorm(), addPeer() and
			state
	*/
	int integrate( SimTime timeEvent, Condition stateEvent = 0 );

protected:
	// Dimension
	int dimension; // highest index

	// integrations steps
	double stepLength;
	double maxStepLength, minStepLength;

	// integration time
	double time;

	// solver error handling
	int relative;
	double errorLimit;

	/**
		\brief Compute error

		\return single error value

		This function is called to compute a single error value
		from the error_vector. The default implementation returns
		the largest error(). A user can overwrite this function to
		use another error norm.

		\sa error, error_value and integrate()
	*/
	virtual double errorNorm();

	// solver
	int stopped;

	std::vector<double> initial_state;
	std::vector<double> slope_1;
	std::vector<double> slope_2;
	std::vector<double> slope_3;

	int reduce();
	void binSearch();
	virtual void takeAStep(double h);

	// stop condition
	bool stopIntegrate();
	SimTime stopTime;
	Condition stopCond;

	// peers
	ProcessList peers;
};

/** \interface ContinuousObserver

	\author Ralf Gerstenberger

	\brief Observer for Continuous specific events

	\sa Continuous

	\since 1.0
*/
class ContinuousObserver:
	public ProcessObserver
{
public:
	virtual void onCreate( Continuous* sender ) {} ///< Construction

	virtual void onAddPeer( Continuous* sender, Process* peer ) {} ///< Add peer
	virtual void onRemovePeer( Continuous* sender, Process* peer ) {} ///< Remove peer

	virtual void onBeginIntegrate( Continuous* sender ) {} ///< Begin integrate
	virtual void onEndIntegrate( Continuous* sender, int result ) {} ///< End integrate

	virtual void onNewValidState( Continuous* sender ) {} ///< New valid state
};

/** \class ContuTrace

	\author Ralf Gerstenberger

	\brief Trace for Continuous state changes

	\sa Continuous ContinuousObserver

	ContuTrace logs the state changes of a Continuous object
	into a text file. The text file is automatically opened
	and closed.

	\since 1.0
*/
class ContuTrace
:	public ContinuousObserver
{
public:
	/// Construction
	ContuTrace( Continuous* contu = 0, const std::string& fileName = "" );
	/// Destruction
	virtual ~ContuTrace();

	/// Follow state changes
	virtual void onNewValidState( Continuous* sender );

private:
	// Implementation
	std::ostream* out;
	bool firstTime;

	void openFile( const std::string& fileName );
};

} } // namespace odemx::base

#endif /* ODEMX_USE_CONTINUOUS */

#endif /* ODEMX_BASE_CONTINUOUS_INCLUDED */
