func.h

a useful spiking neural networks simulator

				}
	virtual real get2ndDeriv(ntime_t t) const
				{ return r_inv_lsq_osq * (
							omega * omega * (- lambda * sin( theta + omega * t ) + omega * cos( theta + omega * t ))
								+ lambda * lambda * exp( lambda * (zerop - t) ) * zerop_value );
				}
	virtual real get3rdDeriv(ntime_t t) const
				{ return r_inv_lsq_osq * (
							omega * omega * omega * (- lambda * cos( theta + omega * t ) - omega * sin( theta + omega * t ))
								- lambda * lambda * lambda * exp( lambda * (zerop - t) ) * zerop_value );
				}

	virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
				{ upslope = fabs(omega); downslope = -upslope; ceil = r; floor = -r;
				  if( getValue(t) > 0 ) {   upslope = get1stDeriv(t); if(   upslope < 0 ) upslope *= -1; }
				                   else { downslope = get1stDeriv(t); if( downslope > 0 ) downslope *= -1; } }

	virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
				{ upslope = omega * omega; downslope = -upslope; ceil = fabs(omega) * r; floor = -ceil;
				  if( get1stDeriv(t) > 0 ) {   upslope = get2ndDeriv(t); if(   upslope < 0 ) upslope *= -1; }
				                      else { downslope = get2ndDeriv(t); if( downslope > 0 ) downslope *= -1; } }

	virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
				{ upslope = fabs( omega * omega * omega); downslope = -upslope; ceil = omega * omega * r; floor = -ceil;
				  if( get2ndDeriv(t) > 0 ) {   upslope = get3rdDeriv(t); if(   upslope < 0 ) upslope *= -1; }
				                      else { downslope = get3rdDeriv(t); if( downslope > 0 ) downslope *= -1; } }

	virtual func_sine_int * clone() { return new func_sine_int(*this); }
	virtual std::string getDescription() { std::ostringstream oss; oss << "func_sine_int r=" << r << ", omega=" << omega << ", theta=" << theta << ", lambda=" << lambda << ", zerop=" << zerop; return oss.str(); }
};

/////////////////////////////////////////////////////////////////////////
/// Integrated function of func_sineshot
///
/// This function represents an integration of func_sineshot.
/// t0 is always equal to zerop.

class func_sineshot_int : public func_deriveq_base
{
	real r;
	real omega;

	real r_div_2;		///< = r/2
	real omega_2;		///< = 2 &omega;
	real duration;		///< = &pi; / &omega;

	static const real alpha = 1.311;
	static const real beta  = 0.375867;

	void valueChange() { duration = M_PI / omega; omega_2 = 2 * omega; r_div_2 = r * 0.5; }
	void zeropChange() {  }

public:
					/// Constructs a func_sineshot_int with amplitude r and angle velocity omega.
	func_sineshot_int(real ir, real iomega) : func_deriveq_base(), r(ir), omega(iomega) { valueChange(); }			

	virtual real getMaxGradient( ntime_t /*t*/ ) const { return fabs(r_div_2) * omega_2; };

	virtual real getValue( ntime_t t ) const
				{ return (zerop <= t && t <= zerop + duration ? r_div_2 * (1-cos(omega_2 * (t-zerop) )) : 0.0); }

	virtual real get1stDeriv( ntime_t t ) const
				{ return (zerop <= t && t <= zerop + duration ? r_div_2 * omega_2 * sin(omega_2 * (t-zerop) ) : 0.0); }

	virtual real get2ndDeriv( ntime_t t ) const
				{ return (zerop <= t && t <= zerop + duration ? r_div_2 * omega_2 * omega_2 * cos(omega_2 * (t-zerop) ) : 0.0); }

	virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
	virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
	virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;

	virtual func_sineshot_int * clone() { return new func_sineshot_int(*this); }
	virtual std::string getDescription() { std::ostringstream oss; oss << "func_sineshot_int r=" << r << ", omega=" << omega; return oss.str(); }
};

/////////////////////////////////////////////////////////////////////////
/// Exponential function
///
/// Exponentially decaying function f(<var>t</var>) = r exp( - &psi; (t - t<sub>0</sub>) ).

class func_exp : public func_base
{
	real r;
	real psi;
	real t0;

public:
			/// Constructs a func_exp with initial value r at time t0, and decay constant psi.
	func_exp(real ir, real ipsi, real it0) : func_base(), r(ir), psi(ipsi), t0(it0) { if( psi < 0.0 ) throw "func_exp"; }	

	virtual real getMaxGradient( ntime_t t ) const { return r>0 ? 0 : get1stDeriv(t); }

	virtual real getValue( ntime_t t ) const
				{ return r * exp( -psi * (t - t0) ); }

	virtual real get1stDeriv( ntime_t t ) const
				{ return -psi * r * exp( -psi * (t - t0) ); }

	virtual real get2ndDeriv( ntime_t t ) const
				{ return psi * psi * r * exp( -psi * (t - t0) ); }

	virtual real get3rdDeriv( ntime_t t ) const
				{ return - psi * psi * psi * r * exp( -psi * (t - t0) ); }

	virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; };

	virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
				{ if( r >= 0 ) {   upslope = 0; ceil =  mak::Infinity; downslope = func_exp::get1stDeriv(t); floor = 0.0; } 
				  else         { downslope = 0; floor= -mak::Infinity;   upslope = func_exp::get1stDeriv(t); ceil  = 0.0; } }

	virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
				{ if( r <= 0 ) {   upslope = 0; ceil =  mak::Infinity; downslope = func_exp::get2ndDeriv(t); floor = 0.0; } 
				  else         { downslope = 0; floor= -mak::Infinity;   upslope = func_exp::get2ndDeriv(t); ceil  = 0.0; } }

	virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
				{ if( r >= 0 ) {   upslope = 0; ceil =  mak::Infinity; downslope = func_exp::get3rdDeriv(t); floor = 0.0; } 
				  else         { downslope = 0; floor= -mak::Infinity;   upslope = func_exp::get3rdDeriv(t); ceil  = 0.0; } }

	virtual func_exp * clone() { return new func_exp(*this); }
	virtual std::string getDescription() { std::ostringstream oss; oss << "func_exp r=" << r << ", psi=" << psi << ", t0=" << t0; return oss.str(); }
};



/////////////////////////////////////////////////////////////////////////
/// Integration of an exponential function
///
/// leaky integration of external, exponential input f(<var>t</var>) = r exp( - &psi; (t - t<sub>0</sub>) ).
/// \todo implement get1stDerivDomain etc. for this function

class func_exp_int : public func_deriveq_base
{
	real r;
	real psi;
	real t0;
	real r_per_lambda_minus_psi;

	void valueChange() { if( psi < 0.0 || lambda < 0.0 ) abort(); r_per_lambda_minus_psi = r / ( lambda - psi ); }

public:
			/// Constructs a func_exp with initial value r at time t0, and decay constant psi.
	func_exp_int(real ir, real ipsi, real it0) : func_deriveq_base(), r(ir), psi(ipsi), t0(it0) { valueChange(); }

	virtual real getMaxGradient( ntime_t /*t*/ ) const { throw "not supported"; }

	virtual real getValue(ntime_t t) const
				{ return r_per_lambda_minus_psi * ( 
								exp( - psi * (t - t0) )
								- exp( lambda * (zerop - t) - psi * (zerop - t0) ) );
				}

	virtual real get1stDeriv(ntime_t t) const
				{ return r_per_lambda_minus_psi * ( 
								- psi * exp( - psi * (t - t0) )
								+ lambda * exp( lambda * (zerop - t) - psi * (zerop - t0) ) );
				}

	virtual real get2ndDeriv(ntime_t t) const
				{ return r_per_lambda_minus_psi * ( 
								psi*psi * exp( - psi * (t - t0) )
								- lambda*lambda * exp( lambda * (zerop - t) - psi * (zerop - t0) ) );
				}

	virtual real get3rdDeriv(ntime_t t) const
				{ return r_per_lambda_minus_psi * ( 
								-psi*psi*psi * exp( - psi * (t - t0) )
								+lambda*lambda*lambda * exp( lambda * (zerop - t) - psi * (zerop - t0) ) );
				}

	virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; };

	virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
				{ 
					real value = getValue(t);
					real deriv = get1stDeriv(t);
					if( r >= 0 )
					{
						if( deriv > 0 )
						{
							real log_lambda_per_psi = log(lambda / psi) / (lambda - psi);
							upslope = deriv; 
							ceil = r_per_lambda_minus_psi * exp( - psi * (zerop - t0) ) * ( exp( -psi * log_lambda_per_psi ) - exp( - lambda * log_lambda_per_psi) );
							downslope = - value / log_lambda_per_psi; /* can be more efficient */
							floor = 0;
						}
						else {
							upslope = 0; ceil = mak::Infinity;
							downslope = -1; throw "Not supported";
							floor = 0;
						}
					}
					else /* r < 0 */ {
						if( deriv < 0 )
						{
							real log_lambda_per_psi = log(lambda / psi) / (lambda - psi);
							downslope = deriv; 
							floor = r_per_lambda_minus_psi * exp( - psi * (zerop - t0) ) * ( exp( -psi * log_lambda_per_psi ) - exp( - lambda * log_lambda_per_psi) );
							upslope = - value / log_lambda_per_psi; /* can be more efficient */
							ceil = 0;
						}
						else {
							downslope = 0; floor = -mak::Infinity;
							upslope = 1; throw "Not supported";
							ceil = 0;
						}
					}
				}

//	virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
	virtual void get1stDerivDomain( ntime_t , real &, real &, real &, real & ) const 
				{ throw "not supported"; }

//	virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
	virtual void get2ndDerivDomain( ntime_t , real &, real &, real &, real & ) const 
				{ throw "not supported"; }

	virtual func_exp_int * clone() { return new func_exp_int(*this); }
	virtual std::string getDescription() { std::ostringstream oss; oss << "func_exp_int r=" << r << ",t0=" << t0; return oss.str(); }
};

/////////////////////////////////////////////////////////////////////////
/// Difference of two exponential functions
///
/// External input of a form f(<var>t</var>) = r1 exp( - &psi;1 (t - t<sub>0</sub>) ) - r2 exp( - &psi;2 (t - t<sub>0</sub>) ).
/// Here r1丄r2 > 0 and 0 < psi1 < psi2.

class func_exp_diff : public func_base
{
	real psi1, psi2;
	real phi, scale, invscale;
	real t0orig;
	real r0, t0;
	real r1, t1;
	real r2, t2;
	real it0, ir1, ir2;
	real r1orig, r2orig;
	real maxval;
	real xi;

	void valueChange();

public:
			/// constructs func_exp_diff.
	func_exp_diff(real ir1, real ipsi1, real ir2, real ipsi2, real it0);

	virtual real getMaxGradient( ntime_t t ) const { return r0<0 ? r0 * xi : (t >= t1 ? 0 : get1stDeriv(t)); }

	virtual real getValue( ntime_t t ) const
				{ return r0 * (exp( -psi1 * (t - t0) ) - exp( -psi2 * (t - t0) ) ); }

	virtual real get1stDeriv( ntime_t t ) const
				{ return r1 * (exp( -psi1 * (t - t1) ) - exp( -psi2 * (t - t1) ) ); }

	virtual real get2ndDeriv( ntime_t t ) const
				{ return r2 * (exp( -psi1 * (t - t2) ) - exp( -psi2 * (t - t2) ) ); }

	virtual real get3rdDeriv( ntime_t t ) const
				{ return r2 * (-psi1 * exp( -psi1 * (t - t2) ) + psi2 * exp( -psi2 * (t - t2) ) ); }

	virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; }

	virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const ;
	virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const ;
	virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const ;

	virtual func_exp_diff * clone() { return new func_exp_diff(*this); }
	virtual std::string getDescription() { std::ostringstream oss; oss << "func_exp_diff r0=" << r0 << ", psi1=" << psi1 << ", psi2=" << psi2 << ", t0=" << t0; return oss.str(); }

	/////////////////////////////////////////////////////////////////////////
	/// Message to add an event. ????

	class message_add_event_time : public message_base
	{
	public:
		message_add_event_time() { }
		static const char * const messageId;
		virtual const char *getMessageId() const { return messageId; }
	};	

	virtual bool processMessage(ntime_t t, const message_base &m);
};


} // namespace punnets

#endif // __func_h