dneuron.cpp

a useful spiking neural networks simulator

				d1d->value = v;
				deriv1st += v;
				(*i)->get1stDerivDomain(current_time, d1d->upslope, d1d->ceil, d1d->downslope, d1d->floor);
				d1d++;
			}
		}

		real d1st_bound = 0.0, d2nd_bound = 0.0;

		if( iporder>=1 && deriv1st >= 1e-5 )
			d1st_bound = domain_floor_time( d1st_domain, next_incontinuity - current_time);
		else if( iporder>=1 && deriv1st <= -1e-5 )
			d1st_bound = domain_ceil_time( d1st_domain, next_incontinuity - current_time);

		if( iporder >= 2 && value_bound < 1e-1 && d1st_bound < 1e-1 && (use_next_known_pulse ? deriv1st + current_time < next_known_pulse : true) )
		{
			vector<valdomain>::iterator d2d = d2nd_domain.begin();
			for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
			{
				real v = (*i)->get2ndDeriv(current_time);
				d2d->value = v;
				deriv2nd += v;
				(*i)->get2ndDerivDomain(current_time, d2d->upslope, d2d->ceil, d2d->downslope, d2d->floor);
				d2d++;
			}

			if( iporder>=2 && deriv2nd >= 1e-5 )
				d2nd_bound = domain_floor_time( d2nd_domain, next_incontinuity - current_time);
			else if( iporder>=2 && deriv2nd <= -1e-5 )
				d2nd_bound = domain_ceil_time( d2nd_domain, next_incontinuity - current_time);
		}
		
		totalpartition++;
		next_bound = 0.0;
		// real bound_val; // The value of the signal at the next bound
		bool need_newton = false;
		if( (iporder<1 || value_bound >= d1st_bound) && (iporder<2 || value_bound >= d2nd_bound) )
		{
			prof<3> profobj("tneuron_ext::3-0th_bound");
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
					<< " Stay, signal level " << signal << ", value_bound is the longest " << value_bound << ":" << d1st_bound << ":" << d2nd_bound << endl; 
			next_bound = value_bound;
			last_simulate_type = 0;
		}
		else if( iporder<2 || d1st_bound >= d2nd_bound )
		{	
			prof<3> profobj("tneuron_ext::3-1st_bound");
			next_bound = d1st_bound;
			last_simulate_type = 1;
			if( deriv1st > 0.0 && (bound_val=calcSignal( current_time + next_bound )) > 0.0 )
				need_newton = true;
//			cout << "signal = " << signal << endl;
//			cout << "deriv1st = " << deriv1st << endl;
//			cout << "bound_val = " << bound_val << endl;
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
					<< " Stay, signal level " << signal << ", d1st_bound is the longest " << value_bound << ":" << d1st_bound << ":" << d2nd_bound << " need_newton = " << (need_newton ? "true" : "false") << endl; 
		}
		else {
			prof<3> profobj("tneuron_ext::3-2nd-order");
			next_bound = d2nd_bound;
			last_simulate_type = 2;
			if( deriv1st < 0.0 && deriv2nd < 0.0 )
			{
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " Stay, signal level " << signal << ", d2nd_bound is the longest " << value_bound << ":" << d1st_bound << ":" << d2nd_bound << ", need_newton = false (d1<0,d2<0)" << endl; 
				need_newton = false;
			}
			else if( (bound_val = calcSignal( current_time + next_bound )) > 0.0 )
			{
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " d2nd_bound is the longest " << value_bound << ":" << d1st_bound << ":" << d2nd_bound << ", need_newton = true (bound_val=" << bound_val << ")" << endl; 
				need_newton = true;
			}
			else if( deriv2nd >= 0.0 )
			{
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " d2nd_bound is the longest " << value_bound << ":" << d1st_bound << ":" << d2nd_bound << ", need_newton = false (no convex nor end)" << endl; 
				need_newton = false;
			}
			else {
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " d2nd_bound is the longest " << value_bound << ":" << d1st_bound << ":" << d2nd_bound << ", peak search" << endl; 
				real left = 0.0, right = next_bound, mid = d1st_bound + (right-d1st_bound) * 0.38197;
				real leftval = signal, rightval = bound_val, midval = calcSignal( current_time + mid );
				if( midval < rightval || midval < leftval )
				{					// Cannot enclose the peak...
					if( midval < rightval ) {
						real right_deriv = 0.0;
						for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
							right_deriv += (*i)->get1stDeriv(current_time + right);
						if( right_deriv >= 0.0 )
						{
							if( getDeb() )
								cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
									<< " peak is not exist, right_deriv = " << right_deriv << endl;
							totalpeaksearch[2]++;
							goto PEAK_SEARCH_FINISH;
						}
					}
					totalpeakenclosing++;
					while( midval < rightval || midval < leftval )
					{
						if( getDeb() )
							cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
								<< " peak enclosing " << setprecision(5) << left << ":" << mid << ":" << right << "=" << leftval << ":" << midval << ":" << rightval << endl;
						if( midval < rightval ) {
							left = mid; leftval = midval; mid = left + (right - left) * 0.38197;
						} else {
							right = mid; rightval = midval; mid = left + (right - left) * 0.61803;
						}
						midval = calcSignal( current_time + mid );
					}
				}
				if( midval > 0.0 )
				{
					if( getDeb() )
						cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
							<< " peak found " << mid << "=" << midval << endl; 
					next_bound = mid; bound_val = midval; need_newton = true;
				}
				else /* Look for the peak */
				{
					real oldp = 0, newp=-1;
					while( fabs(oldp - newp) >= 1e-8 ) 
					{
						if( getDeb() )
							cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
								<< " peak search " << setprecision(5) << left << ":" << mid << ":" << right << "=" << leftval << ":" << midval << ":" << rightval << endl;

						real newpmove, newpval;
						real mlmr = (mid - left) * (midval - rightval);
						real mrml = (mid - right) * (midval - leftval);

						// abort peak search, if no hope to go over zero

						real left_best = midval - mlmr / (mid - right);
						real right_best = midval - mrml / (mid - left);
						if( left_best < 0 && right_best < 0 )
						{
							if( getDeb() )
								cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
									<< " peak search aborted, left_best=" << setprecision(10) << left_best << ", right_best=" << right_best << endl;
							break;
						}

						real mlmr_mrml = -2.0 * (mlmr - mrml);
						real mlmlmr_mrmrml = (mid - left) * mlmr - (mid - right) * mrml;
						if( mlmr_mrml < 0.0 ) { mlmr_mrml = - mlmr_mrml; mlmlmr_mrmrml = - mlmlmr_mrmrml; }
						
						if( mlmlmr_mrmrml > mlmr_mrml*(left-mid) && mlmlmr_mrmrml < mlmr_mrml*(right-mid) )
						{
							newpmove = mlmlmr_mrmrml / mlmr_mrml; // use parabora interpolation
							if( getDeb() )
								cout << "parabola " << setprecision(10) << newpmove << ", " << newp - (mid + newpmove) << endl;
						}
						else
						{
							cout << "mid - left = " << setprecision(10) << mid - left << endl;
							cout << "mid - right = " << setprecision(10) << mid - right << endl;
							cout << "parabola = " << setprecision(10) << mlmlmr_mrmrml / mlmr_mrml << endl;
							if( mid - left > right - mid ) newpmove = (left-mid) * 0.38197;
													  else newpmove = (right-mid) * 0.38197;
							cout << "golden   " << setprecision(10) << newpmove << ", " << newp - (mid + newpmove) << endl;
						}
						oldp = newp;
						newp = mid + newpmove;
						newpval = calcSignal(current_time + newp);
						if( newpval > midval )
						{ 
							if( newp < mid ) { right = mid; rightval = midval; }
							            else { left = mid; leftval = midval; }
							mid = newp; midval = newpval; 
							if( newpval > 0.0 )
							{ next_bound = newp; bound_val = newpval; need_newton = true; break; }
						}
						else 
							if( newp < mid ) { left = newp; leftval = newpval; }
							            else { right = newp; rightval = newpval; }
					}
					/* if peak value is < 0.0, then need_newton is false */
					totalpeaksearch[need_newton ? 1 : 0]++;
				}
				PEAK_SEARCH_FINISH:
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " peak search finished, need_newton = " << (need_newton ? "true" : "false") << ", mid " << mid << "=" << midval << endl;
			}
		}
		if( need_newton == false )
		{
			if( next_bound < ip_delta_t )
			{
				next_bound = ip_delta_t;
				last_simulate_type = 3;
				if( (bound_val=calcSignal( current_time + next_bound )) > 0.0 )
					need_newton = true;
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " next_bound is adjusted to delta_t, need_newton = " << (need_newton ? "true" : "false") << endl;
			}
			if( need_newton == false ) {
				totalpartition_nonewton[last_simulate_type]++;
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " Reschedule at the end of no-newton-area " << current_time + next_bound << endl; 
				setLoopBack(scheduler, current_time + next_bound);
				return; 
			}
		}
		}
		prof<2> profobj2("tneuron_ext::4 newton");
		totalpartition_newton[last_simulate_type]++;
		/* Do newton to find fire point */
		if( next_bound >= 10000000 )
		{
			next_bound = 10000000; // to avoid infinity operation
			if( bound_val >= 10000000 )
				bound_val = calcSignal(current_time + next_bound);
		}
		real left = 0.0, right = next_bound;
		real newtonp = right * signal / (signal - bound_val); // set newtonp to the linear interpolation point
		if( getDeb() )
			cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
				<< " Start newton in range of " << next_bound << ", signal=" << signal << ", bound_val=" << bound_val << endl; 

		real val = 0.0, deriv= 0.0;
		for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
		{
			val += (*i)->getValue(current_time + newtonp);
			deriv += (*i)->get1stDeriv(current_time + newtonp);
		}

		real old_newtonp = newtonp - 1;
		while( fabs(old_newtonp - newtonp) > 1e-8 )
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
					<< " Newton left=" << left << ", newtonp=" << newtonp << ", right=" << right << ", val=" << val << endl;
			old_newtonp = newtonp;

			real newtonp_move = Infinity;
			if( fabs(deriv) > 1e-8 ) 
				newtonp_move = - val / deriv;
			if( newtonp + newtonp_move < left || newtonp + newtonp_move > right )
			{
				if( getDeb() )
					cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
						<< " Newton_safe newtonp_move=" << newtonp + newtonp_move << endl;
				if( val > 0.0 ) newtonp_move = (left - newtonp) * 0.5;
				           else newtonp_move = (right - newtonp) * 0.5;
			}
			if( newtonp_move < 0.0 )
				{ right = newtonp; newtonp += newtonp_move; }
			else{ left  = newtonp; newtonp += newtonp_move; }

			val = 0.0; deriv= 0.0;
			for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
			{
				val += (*i)->getValue(current_time + newtonp);
				deriv += (*i)->get1stDeriv(current_time + newtonp);
			}
		}
		if( val > -1e-5 )
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
					<< " Newton success, reschedule at next fire " << current_time + newtonp << endl; 
			setLoopBack(scheduler, current_time + newtonp);
		}
		else
		{
			if( getDeb() )
				cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
					<< " Newton failed, Reschedule at next prediction bound " << current_time + next_bound << endl; 
			setLoopBack(scheduler, current_time + next_bound);
		}
	}
//	if( getDeb() ) cout << "leave tneuron_ext::schedulefire" << endl;
}

// template instantiation.
template class punnets_private::tsynapse<true>        ;
template class punnets_private::tsynapse_message<true>;
template class punnets_private::tsynapse_fatigue<true>;
template class punnets_private::tsynapse_addfunc<true>   ;
template class punnets_private::tsynapse_messfunc<true>   ;

template class punnets_private::tneuron<true>          ;
template class punnets_private::tneuron_ext_const<true>;
template class punnets_private::tneuron_ext<true>      ;

template class punnets_private::tsynapse<false>        ;
template class punnets_private::tsynapse_message<false>;
template class punnets_private::tsynapse_fatigue<false>;
template class punnets_private::tsynapse_addfunc<false>   ;
template class punnets_private::tsynapse_messfunc<false>   ;

template class punnets_private::tneuron<false>          ;
template class punnets_private::tneuron_ext_const<false>;
template class punnets_private::tneuron_ext<false>      ;

/////////////////////////////////////////////////////////////////////////
} // namespace punnets_private
/////////////////////////////////////////////////////////////////////////