func.cpp

a useful spiking neural networks simulator

/** 
 * @file  func.cpp
 * @brief Functions
 *
 * @author Makino, Takaki <[e-mail]>
 * @date 2003-05-01
 * @version $Id: func.cpp,v 1.3 2003/05/02 09:14:21 t Exp $
 *
 *  Copyright (C) 2003 Makino, Takaki. 
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2, or (at your option)
 *  any later version.

 *  This program 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.  
 */

#include "func.h"
#include <iostream>

namespace punnets_common {

using namespace mak;


const char * const func_deriveq_base::message_set_lambda::messageId = "message_set_lambda";
const char * const func_deriveq_base::message_set_zero_point::messageId = "message_set_zerop";
const char * const func_delta_int::message_add_pulse::messageId = "message_add_pulse";
const char * const func_sineshot::message_set_t0::messageId = "message_set_t0";
const char * const func_exp_diff::message_add_event_time::messageId = "message_add_event_time";

/***************************************************************************/
/**************** func_const_int *******************************************/
/***************************************************************************/

void func_const_int::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	if( c > 0 ) { upslope = func_const_int::get1stDeriv(t); downslope = 0; ceil = limit; floor = 0;  }
	else        { downslope = func_const_int::get1stDeriv(t); upslope = 0; floor = limit; ceil = 0;  }
}

void func_const_int::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	if( c > 0 ) { upslope = 0; downslope = func_const_int::get2ndDeriv(t); ceil = limit * lambda; floor = 0;  }
	else        { downslope = 0; upslope = func_const_int::get2ndDeriv(t); floor = limit * lambda; ceil = 0;  }
}

void func_const_int::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	if( c > 0 ) { upslope = func_const_int::get3rdDeriv(t); downslope = 0; ceil = 0; floor = -Infinity;  }
	else        { downslope = func_const_int::get3rdDeriv(t); upslope = 0; floor = 0; ceil = Infinity;  }
}

/***************************************************************************/
/**************** func_delta_int *******************************************/
/***************************************************************************/

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

void func_delta_int::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	if( zerop >= t0 || t0 > t ) { upslope = downslope = 0; ceil = floor = 0; }  
	else if( r <= 0 ) {   upslope = 0; ceil =  Infinity; downslope = get2ndDeriv(t); floor = 0.0; } 
	else              { downslope = 0; floor= -Infinity;   upslope = get2ndDeriv(t); ceil  = 0.0; } 
}

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

/***************************************************************************/
/**************** func_response *******************************************/
/***************************************************************************/

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

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

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

/***************************************************************************/
/**************** func_sine ************************************************/
/***************************************************************************/

void func_sine::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	real phi = omega * t + theta;
	if( phi < 0.0 ) phi -= ::floor(phi / (2*M_PI)) * 2*M_PI;
	phi = fmod(phi, 2*M_PI);
	phi -= M_PI;
	if( phi + 0.5*M_PI >= 0 )
		upslope = omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
	else
		upslope = omega * r * cos(omega * t + theta);
	phi += M_PI;
	if( phi >= M_PI ) phi -= 2*M_PI;
	if( phi + 0.5*M_PI >= 0 )
		downslope = - omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
	else
		downslope = omega * r * cos(omega * t + theta);
	ceil = r; floor = -r; 
}

void func_sine::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	real phi = omega * t + theta + 0.5*M_PI;
	if( phi < 0.0 ) phi -= ::floor(phi / (2*M_PI)) * 2*M_PI;
	phi = fmod(phi, 2*M_PI);
	phi -= M_PI;
	if( phi + 0.5*M_PI >= 0 )
		upslope = omega * omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
	else
		upslope = -omega * omega * r * cos(phi);
	phi += M_PI;
	if( phi >= M_PI ) phi -= 2*M_PI;
	if( phi + 0.5*M_PI >= 0 )
		downslope = - omega * omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
	else
		downslope = omega * omega * r * cos(phi);
	ceil = r*omega; floor = -r*omega; 
}

void func_sine::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{
	real phi = omega * t + theta + M_PI;
	if( phi < 0.0 ) phi -= ::floor(phi / (2*M_PI)) * 2*M_PI;
	phi = fmod(phi, 2*M_PI);
	phi -= M_PI;
	if( phi + 0.5*M_PI >= 0 )
		upslope = omega * omega * omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
	else
	    upslope = -omega * omega * omega * r * cos(phi);
	phi += M_PI;
	if( phi >= M_PI ) phi -= 2*M_PI;
	if( phi + 0.5*M_PI >= 0 )
		downslope = - omega * omega * omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
	else
		downslope = omega * omega * omega * r * cos(phi);
	ceil = r*omega*omega; floor = -r*omega*omega; 
}


/***************************************************************************/
/**************** func_sineshot ********************************************/
/***************************************************************************/

void func_sineshot::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	upslope = ceil = downslope = floor = 0.0;
	if( t < t0 || t >= t0 + duration )
		return;
	else 
	{
		real phase = omega_2 * (t - t0);
		real phi = phase + M_PI;
		if( phi >= 2*M_PI ) phi -= 2*M_PI;
		( r > 0 ? ceil : floor ) = r;
		if( phase < 0.5 * M_PI )
		{
			( r > 0 ? upslope : downslope ) = omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
			( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi-M_PI) ) + 0.5*M_PI - alpha);
		}
		else if( phase < M_PI )
		{
			( r > 0 ? upslope : downslope ) = r_div_2 * omega_2 * sin(phase );
			( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi-M_PI) ) + 0.5*M_PI - alpha);
		}
		else if( phase < 1.5*M_PI )
		{
			( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi+M_PI) ) + 0.5*M_PI - alpha);
		}
		else
		{
			( r > 0 ? downslope : upslope ) = r_div_2 * omega_2 * sin(phase);
		}
	}
}

void func_sineshot::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	upslope = ceil = downslope = floor = 0.0;
	if( t < t0 || t >= t0 + duration )
		return;
	else 
	{
		real phase = omega_2 * (t - t0);
		real phi = phase + 1.5*M_PI;
		if( phi >= 2*M_PI ) phi -= 2*M_PI;
		ceil = fabs(r_div_2) * omega_2;
		floor = -ceil;
		if( phase < 0.5 * M_PI )
		{
			( r > 0 ? upslope : downslope ) = r_div_2 * omega_2 * omega_2 * cos(phase );
			( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi - M_PI) ) + 0.5*M_PI - alpha);
		}
		else if( phase < M_PI )
		{
		//	( r > 0 ? upslope : downslope ) = 0; 
			( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi + M_PI) ) + 0.5*M_PI - alpha);
		}
		else if( phase < 1.5*M_PI )
		{
			( r > 0 ? ceil : floor ) = 0; 
			( r > 0 ? upslope : downslope ) = (- r_div_2 * omega_2 * sin(phase )) / (duration - (t-t0));
			( r > 0 ? downslope : upslope ) = r_div_2 * omega_2 * omega_2 * cos(phase );
		}
		else
		{
			( r > 0 ? ceil : floor ) = 0; 
			( r > 0 ? upslope : downslope ) = (- r_div_2 * omega_2 * sin(phase )) / (duration - (t-t0));
		}
	}
}

void func_sineshot::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	upslope = ceil = downslope = floor = 0.0;
	if( t < t0 || t >= t0 + duration )
		return;
	else 
	{
		real phase = omega_2 * (t - t0);
		real phi = phase;
		ceil = fabs(r_div_2 * omega_2 * omega_2);
		floor = -ceil;
		if( phase < 0.5 * M_PI )
		{
			( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
			( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi + M_PI) ) + 0.5*M_PI - alpha);
		}
		else if( phase < M_PI )
		{
			( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
			( r > 0 ? downslope : upslope ) = - r_div_2 * omega_2 * omega_2 * omega_2 * sin( phase );
 		}
		else if( phase < 1.5*M_PI )
		{
			( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
		}
		else
		{
			( r > 0 ? upslope : downslope ) = - r_div_2 * omega_2 * omega_2 * omega_2 * sin( phase );
			( r > 0 ? downslope : upslope ) = - (r_div_2 * omega_2 * omega_2 * cos(phase )) / (duration - (t-t0));
			( r > 0 ? floor : ceil ) = 0; 
		}
	}