func.cpp

a useful spiking neural networks simulator

}


/***************************************************************************/
/**************** func_sineshot_int ********************************************/
/***************************************************************************/

void func_sineshot_int::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	const real t0 = zerop;
	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_int::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	const real t0 = zerop;
	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_int::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	const real t0 = zerop;
	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; 
		}
	}
}

/***************************************************************************/
/**************** func_exp_diff ********************************************/
/***************************************************************************/

func_exp_diff::func_exp_diff(real iir1, real ipsi1, real iir2, real ipsi2, real iit0) 
	: func_base(), psi1(ipsi1), psi2(ipsi2), t0orig(iit0), r1orig(iir1), r2orig(iir2)
{
	it0 = t0orig;
	ir1 = r1orig;
	ir2 = r2orig;
	valueChange();
};

void func_exp_diff::valueChange()
{
	if( psi1 < 0.0 || psi2 <= psi1 || ir1 *ir2 <= 0.0 ) 
	{
		std::cout << "Abort!! bad condition!! " << std::endl;
		std::cout << "psi1 = " << psi1 << ", psi2 = "<< psi2 << ", ir1 = " << ir1 << ", ir2 = " << ir2 << std::endl;
		throw "func_exp_diff"; 
	}

	phi = psi2 / psi1;
	real z = 2 * log( phi ) / (phi - 1);
	scale = psi1 / z;
	invscale = z / psi1;
	maxval = exp(-psi1 * 0.5 * invscale) - exp(-psi2 * 0.5 * invscale);
	xi = - psi1 * exp(-psi1 * invscale) + psi2 * exp(-psi2 * invscale);

	t0 = (log(ir1/ir2))/(psi1-psi2);
	r0 = ir1 * exp( -psi1 * t0 ); 

	t1 = (log(ir1*psi1/(ir2*psi2)))/(psi1-psi2);
	r1 = - ir1 * psi1 * exp( -psi1 * t1 );

	t2 = (log(ir1*psi1*psi1/(ir2*psi2*psi2)))/(psi1-psi2);
	r2 = ir1 * psi1 * psi1 * exp( -psi1 * t2 );

	if( fabs( ir1 * exp( -psi1 * t0 ) - ir2 * exp( -psi2 * t0 ) ) > 1e-5 ) 
		std::cout << "komatta 0! " << ir1 * exp( -psi1 * t0 ) << ", " << ir2 * exp( -psi2 * t0 ) << std::endl << getDescription() << std::endl;
	if( fabs( ir1 * psi1 * exp( -psi1 * t1 ) - ir2 * psi2 * exp( -psi2 * t1 ) ) > 1e-5 )
		std::cout << "komatta 1! " << ir1 * psi1 * exp( -psi1 * t1 ) << ", " << ir2 * psi2 * exp( -psi2 * t1 ) << std::endl;
	if( fabs( ir1 * psi1 * psi1 * exp( -psi1 * t2 ) - ir2 * psi2 * psi2 * exp( -psi2 * t2 ) ) > 1e-5 )
		std::cout << "komatta 2! " << ir1 * psi1 * psi1 * exp( -psi1 * t2 ) << ", " << ir2 * psi2 * psi2 * exp( -psi2 * t2 ) << std::endl;
//	cout << "scale = "<< scale << ", maxval = "<< maxval << ", xi= "<< xi << endl;
//	cout << "t0 = "<< t0 << ", r0 = "<< r0 << endl;
//	cout << "t1 = "<< t1 << ", r1 = "<< r1 << endl;
//	cout << "t2 = "<< t2 << ", r2 = "<< r2 << endl;
	t0 += it0; t1 += it0; t2 += it0;
	phi = 0.3 *(log(phi)/log(10)) * (log(phi)/log(10)) + 2.45;
}	

void func_exp_diff::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{
	real val = getValue(t);
	real der = get1stDeriv(t);

	if( t < t0 + 0.5 * invscale ) 
	{
		if( t < t0 ) 
		{
			if( r0 > 0 ) { ceil  = maxval * r0; floor = val; upslope   = der; downslope = 0; }
			else         { floor = maxval * r0; ceil  = val; downslope = der; upslope = 0;   }
		}
		else // if( t < t0 + 0.5 * invscale ) 
		{
			if( r0 > 0 ) { ceil  = maxval * r0; floor = 0; upslope   = der; downslope = r0 * xi * (1.0 - pow(1.0 - (t-t0)*scale, phi)); }
			else         { floor = maxval * r0; ceil  = 0; downslope = der; upslope   = r0 * xi * (1.0 - pow(1.0 - (t-t0)*scale, phi)); }
		}
	}
	else if( t < t0 + 1.0 * invscale ) 
	{
		if( r0 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = r0 * xi * (1.0 - pow(1.0 - (t-t0)*scale, phi)); }
		else         { floor = val; ceil  = 0; downslope = 0; upslope   = r0 * xi * (1.0 - pow(1.0 - (t-t0)*scale, phi)); }
	}
	else
	{
		if( r0 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = der; }
		else         { floor = val; ceil  = 0; downslope = 0; upslope   = der; }
	}
}

void func_exp_diff::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
	real val = get1stDeriv(t);
	real der = get2ndDeriv(t);

	if( t < t1 + 0.5 * invscale ) 
	{
		if( t < t1 ) 
		{
			if( r1 > 0 ) { ceil  = maxval * r1; floor = val; upslope   = der; downslope = 0; }
			else         { floor = maxval * r1; ceil  = val; downslope = der; upslope = 0;   }
		}
		else // if( t < t1 + 0.5 * invscale ) 
		{
			if( r1 > 0 ) { ceil  = maxval * r1; floor = 0; upslope   = der; downslope = r1 * xi * (1.0 - pow(1.0 - (t-t1)*scale, phi)); }
			else         { floor = maxval * r1; ceil  = 0; downslope = der; upslope   = r1 * xi * (1.0 - pow(1.0 - (t-t1)*scale, phi)); }
		}
	}
	else if( t < t1 + 1.0 * invscale ) 
	{
		if( r1 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = r1 * xi * (1.0 - pow(1.0 - (t-t1)*scale, phi)); }
		else         { floor = val; ceil  = 0; downslope = 0; upslope   = r1 * xi * (1.0 - pow(1.0 - (t-t1)*scale, phi));   }
	}
	else
	{
		if( r1 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = der; }
		else         { floor = val; ceil  = 0; downslope = 0; upslope   = der; }
	}
}

void func_exp_diff::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{
	real val = get2ndDeriv(t);
	real der = get3rdDeriv(t);

	if( t < t2 + 0.5 * invscale ) 
	{
		if( t < t2 ) 
		{
			if( r2 > 0 ) { ceil  = maxval * r2; floor = val; upslope   = der; downslope = 0; }
			else         { floor = maxval * r2; ceil  = val; downslope = der; upslope = 0;   }
		}
		else // if( t < t2 + 0.5 * invscale ) 
		{
			if( r2 > 0 ) { ceil  = maxval * r2; floor = 0; upslope   = der; downslope = r2 * xi * (1.0 - pow(1.0 - (t-t2)*scale, phi)); }
			else         { floor = maxval * r2; ceil  = 0; downslope = der; upslope   = r2 * xi * (1.0 - pow(1.0 - (t-t2)*scale, phi)); }
		}
	}
	else if( t < t2 + 1.0 * invscale ) 
	{
		if( r2 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = r2 * xi * (1.0 - pow(1.0 - (t-t2)*scale, phi)); }
		else         { floor = val; ceil  = 0; downslope = 0; upslope   = r2 * xi * (1.0 - pow(1.0 - (t-t2)*scale, phi));   }
	}
	else
	{
		if( r2 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = der; }
		else         { floor = val; ceil  = 0; downslope = 0; upslope   = der; }
	}
}

bool func_exp_diff::processMessage(ntime_t t, const message_base &m) 
{ 
#ifdef USE_DYNAMIC
	if( dynamic_cast<const message_add_event_time *>(&m) != NULL ) 
#else
	if( m.getMessageId() == message_add_event_time::messageId ) 
#endif
	{
		ir1 = r1orig + ir1 * exp( - psi1 * (t-it0) );
		ir2 = r2orig + ir2 * exp( - psi2 * (t-it0) );
		it0 = t;
		valueChange();
		return true;
	}
	return func_base::processMessage(t, m); 
}


} // namespace punnets