---

dneuron.cpp

Go to the documentation of this file.

#include "dneuron.h"
#include <math.h>
#include <mak/profile.h>

using namespace mak;

//#define DEBUG

namespace punnets_common {

u_int64_t totalfire = 0;
u_int64_t totalpulse = 0;
u_int64_t totalpartition = 0;
u_int64_t totalpartition_nonewton[4] = {0,0,0,0};
u_int64_t totalpartition_newton[4] = {0,0,0,0};
u_int64_t totalpeaksearch[3] = {0,0,0};
u_int64_t totalpeakenclosing = 0;
u_int64_t totalrescheduled = {0};
u_int64_t totalfiltered_maxgrad = 0;
u_int64_t totalfiltered_incontinuity = 0;
u_int64_t totalfiltered_nextpulse = 0;

const int iporder = 2;
const ntime_t ip_delta_t = 1e-5;

//const bool use_next_known_pulse = false;
const bool use_next_known_pulse = true;

using namespace std;

class tpulse : public taction
{
    tneuron_base &dest;
    real level;

public:
    tpulse(tneuron_base &idest, real ilevel) : taction(), dest(idest), level(ilevel) { }
    
    virtual void activate(tscheduler &scheduler, ntime_t current_time) 
        { 
//      if( getDeb() )
//          cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << dest.getName() << " Pulse Arrived (external input)" << endl;
        dest.pulseArrive( scheduler, current_time, level ); delete this; }

    tqueue *queue() const { return dest.queue(); }
    virtual const char *getClassName() const { return "tpulse"; };
};

class tmessage : public taction
{
    tneuron_base &dest;
    message_base *mess;
    real level;

public:
    tmessage(tneuron_base &idest, message_base *imess) : taction(), dest(idest), mess(imess) { }
    virtual ~tmessage() { delete mess; }
    
    virtual void activate(tscheduler &scheduler, ntime_t current_time) 
        { 
//      if( getDeb() )
//          cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << dest.getName() << " Pulse Arrived (external message " << mess->getMessageId() << ")" << endl;
        dest.pulseArrive( scheduler, current_time, mess ); delete this; }

    tqueue *queue() const { return dest.queue(); }
    virtual const char *getClassName() const { return "tmessage"; };
};


taction &makePulse( tneuron_base &idest, real ilevel )
{
    return *new tpulse( idest, ilevel );
}

taction &makePulse( tneuron_base &idest, message_base *mess )
{
    return *new tmessage( idest, mess );
}


}  // namespace punnets

namespace punnets_private {

template <bool debug>
tneuron<debug>::~tneuron()

{
    for( vector<tsynapse_base *>::iterator i = synapses.begin(); i != synapses.end(); i++ )
        delete *i;
}


template <bool debug>
tneuron<debug>::tneuron(string iname /*= ""*/, 
                 ntime_t isig_hv_period /*= def_sig_hv_period */,
                 ntime_t ithr_hv_period /*= def_thr_hv_period */,
                 real imin_threshold    /*= def_min_threshold */,
                 real imax_threshold    /*= def_max_threshold */ )
        : tneuron_base(iname), 
          sig_level(0.0), last_simulate(0.0), last_fire(0.0),
          sig_hv_period(isig_hv_period), thr_hv_period(ithr_hv_period),
          min_threshold(imin_threshold), max_threshold(imax_threshold),
          coeff_sigdecay(- M_LN2 / isig_hv_period), 
          coeff_thrdecay( - M_LN2 / ithr_hv_period )
{ }


template <bool debug>
void tneuron<debug>::fire(tscheduler &scheduler, ntime_t current_time)
{
    prof<1> pf("tneuron::Fire");
    totalfire++;

    for( vector<tsynapse_base *>::iterator i = synapses.begin(); i != synapses.end(); i++ )
    {
//#ifdef DEBUG
//      cout << "schedule pulse at " << current_time + i->getDelay() << " to " << i->getDest().getName() << endl;
//#endif
        scheduler.scheduleEvent(current_time + (*i)->getDelay(), **i);
    }
    sig_level = 0.0;
    last_fire = current_time;
}

template <bool debug>
void tneuron<debug>::scheduleFire(tscheduler &scheduler, ntime_t current_time, bool /*resched*/)
{
    prof<1> pf("tneuron::scheduleFire");
    real current_threshold = 0.0;
    if( sig_level >= min_threshold )
    {
        // �ő債‚«‚¢’l‚ð’´‚¦‚Ä‚¢‚½‚瑦�À‚É”­‰Î
        if( sig_level >= max_threshold || (current_threshold = getCurrentThrLevel( current_time )) - sig_level < epsilon )
        {
            if( getDeb() )
                cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Fire, level " << sig_level << ", threshold " << current_threshold << endl;
            fire(scheduler, current_time);
            return;
        }
        else
        {
            if( getDeb() )
                cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Candid, level " << sig_level << ", threshold " << current_threshold << endl;
        }
    }
    else
        if( getDeb() )
            cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Accumulate, " << sig_level << endl;
}

/************************************************************************/

template <bool debug>
void tneuron_ext_const<debug>::scheduleFire(tscheduler &scheduler, ntime_t current_time, bool resched)
{
    real current_threshold = 0.0;
    if( sig_level >= min_threshold || sig_converge_level > min_threshold )
    {
        // �ő債‚«‚¢’l‚ð’´‚¦‚Ä‚¢‚½‚瑦�À‚É”­‰Î
        if( sig_level >= max_threshold || (current_threshold = getCurrentThrLevel( current_time )) - sig_level < epsilon )
        {
            if( getDeb() )
                cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Fire, level " << sig_level << ", threshold " << current_threshold << endl;
            fire(scheduler, current_time);
            if( sig_converge_level < min_threshold )
                return;
            current_threshold = max_threshold;
        }
        else
        {
            if( getDeb() )
                cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Candid, level " << sig_level << ", threshold " << current_threshold << endl;
            if( ! resched )
                return;
        }

        // �¡Œã�A‚µ‚«‚¢’lˆÈ‰º‚Å‚ ‚ê‚Î�A’´‚¦‚é‰Â”\�«‚ª‚ ‚é‚©’²‚ׂé
        real sig_decay_level = sig_level - sig_converge_level;
        real thr_decay_level = current_threshold - min_threshold;
        // �Å�¬‚µ‚«‚¢’l‚ª�M�†Žû‘©’l–¢–ž‚È‚çŠmŽÀ‚ɉz‚¦‚é‚Ì‚Å�A‚»‚êˆÈŠO‚Ì�ê�‡‚¾‚¯’²‚ׂé
        if( min_threshold >= sig_converge_level )
        {
            // �M�†‚Ì”¼Œ¸Šú‚Ì‚Ù‚¤‚ª’·‚¢�ê�‡ˆÈŠO‚͉”\�«‚ª‚È‚¢
            if( coeff_sigdecay - coeff_thrdecay <= 0 )
            {
                if( getDeb() )
                    cout << current_time << " Fail neuron " << name << endl;
                return;
            }
                
            // (‚µ‚«‚¢’lŠÖ�” - �M�†ŠÖ�”) ‚ª‹É‘å‚ɂȂ鎞��
            ntime_t maximal_time = 
                (coeff_sigdecay - coeff_thrdecay) *
                log( ( coeff_thrdecay * thr_decay_level ) / 
                     ( coeff_sigdecay * sig_decay_level )   );
            // ‹É‘å‚É‚È‚Á‚½‚Æ‚«‚ɉz‚¦‚Ä‚¢‚é‚©�H
            real siglvl_at_maximal = sig_decay_level * exp(coeff_sigdecay * maximal_time) + sig_converge_level;
            real thrlvl_at_maximal = thr_decay_level * exp(coeff_thrdecay * maximal_time) + min_threshold;
            if( siglvl_at_maximal < thrlvl_at_maximal )
            {
                if( getDeb() )
                    cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Fail, max_lvl " << siglvl_at_maximal << ", threshold " << thrlvl_at_maximal << endl;
                return; // ‰z‚¦‚Ä‚¢‚È‚¢�A”­‰Î‚̉”\�«‚È‚µ
            }
            if( getDeb() )
                cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Possible, max_lvl " << siglvl_at_maximal << ", threshold " << thrlvl_at_maximal << endl;
//          cout << T << endl;
        }
        else
        {
            if( getDeb() )
                cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Possible" << endl;
//          cout << T << endl;
        }
        // ‚µ‚΂炭Œã‚Å”­‰Î‚·‚é‚Ì‚Å‚»‚ÌŽž�� t ‚𒲂ׂé
        // ƒjƒ…�[ƒgƒ“–@ T = exp(t)

        // �‰Šú’l‚Í 1.0 ‚¾‚ª�A‚P‰ñ‚܂킵‚½Œã‚Ì’l‚ÅŒvŽZ
        real T = 1.0 - (sig_level - current_threshold) / (coeff_sigdecay * sig_decay_level - coeff_thrdecay * thr_decay_level);

        int count = 0;

        ntime_t fire_time = log(T);
        while(1) {
            real atb = sig_decay_level * exp(fire_time * coeff_sigdecay);
            real ctd = thr_decay_level * exp(fire_time * coeff_thrdecay);
            real T_ = (atb - ctd - min_threshold + sig_converge_level) / (coeff_sigdecay * atb - coeff_thrdecay * ctd);
            if( fabs(T_ ) < epsilon )
                break;
            T *= 1.0 - T_;
//              cout << T << "\t" << T_ << endl;
            fire_time = log(T);
            count++;
        }

        // ”­‰ÎŽž��‚É‚à‚¤ˆê“xƒXƒPƒWƒ…�[ƒŠƒ“ƒO
        if( getDeb() )
        {
            real siglvl_at_fire = sig_decay_level * exp(coeff_sigdecay * fire_time) + sig_converge_level;
            real thrlvl_at_fire = thr_decay_level * exp(coeff_thrdecay * fire_time) + min_threshold;
            cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Reschedule, time " << fire_time << ", fire_lvl " << siglvl_at_fire << ", threshold " << thrlvl_at_fire << ", newton " << count << endl;
        }
//              cout << "newton " << count << endl;
        scheduler.scheduleEvent(current_time + fire_time, *this);
    }
    else
        if( getDeb() )
            cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name << " Accumulate, level " << sig_level << endl;
}

/************************************************************************/


/************************************************************************
 ***** tneuron_ext ******************************************************
 ************************************************************************/

template <bool debug>
tneuron_ext<debug>::~tneuron_ext()
{
    for( vector<tsynapse_base *>::iterator i = synapses.begin(); i != synapses.end(); i++ )
        delete *i;
    for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
        delete *i;
}

template <bool debug>
tneuron_ext<debug>::tneuron_ext(string iname /*= ""*/, real sig_hv_period)
        : tneuron_base(iname), 
          last_simulate(-10000000), last_fire(-10000000),
          lambda(M_LN2 / sig_hv_period) 
{
    pulses = new func_delta_int(0, 0);
    addExt(pulses);
}

template <bool debug>
bool tneuron_ext<debug>::sendMessage(ntime_t t, message_base *mess)
{
    for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
        if( (*i)->processMessage(t, *mess) )
            return true;
    return false;
}

template <bool debug>
bool tneuron_ext<debug>::broadcastMessage(ntime_t t, message_base *mess)
{
    bool ret = false;
    for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
        ret = (*i)->processMessage(t, *mess) || ret;
    return ret;
}

template <bool debug>
real tneuron_ext<debug>::calcSignal( ntime_t current_time )
{
    real signal = 0.0;
    for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
    {
        signal += (*i)->getValue(current_time);
    }
    return signal;
}

template <bool debug>
void tneuron_ext<debug>::fire(tscheduler &scheduler, ntime_t current_time)
{
    {
        prof<1> profobj("Fire");
        totalfire++;

        for( vector<tsynapse_base *>::iterator i = synapses.begin(); i != synapses.end(); i++ )
        {
    //#ifdef DEBUG
    //      cout << "schedule pulse at " << current_time + i->getDelay() << " to " << i->getDest().getName() << endl;
    //#endif
            scheduler.scheduleEvent(current_time + (*i)->getDelay(), **i);
        }
        last_fire = current_time;
        for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
        {
            (*i)->setZeroPoint(current_time);
        }
    }
    scheduleFire(scheduler, current_time);
}


struct valdomain { real value, upslope, ceil, downslope, floor; };

ostream& operator<<(ostream &o, const valdomain &v) 
{
    return o << setprecision(8) << "val=" << v.value << ", u=" << v.upslope << ", c=" << v.ceil << ", d=" << v.downslope << ", f=" << v.floor;
}

const bool debug_domain = false;

real domain_floor_time( vector<valdomain> &sig_domain, ntime_t next_incontinuity )
{
    if( debug_domain )
        cout << "enter domain_floor" << endl;
    ntime_t search_from = 0;
    ntime_t old_search_from;

    do {
        if( debug_domain )
            cout << "domain_floor search_from=" << search_from << endl;

        real sig = 0.0, slo = 0.0;
        for( vector<valdomain>::iterator i = sig_domain.begin(); i != sig_domain.end(); i++ )
        {
            real nsig = i->value + i->downslope * search_from;
            if( debug_domain )
                cout << "domain nsig=" << nsig << ", " << *i << endl;
            if( nsig < i->floor ) {
                sig += i->floor;
            }
            else {
                sig += nsig; 
                slo += i->downslope;
            }
        }
        old_search_from = search_from;
        if( debug_domain )
            cout << "domain_floor sig=" << sig << ", slo=" << slo << endl;
        if( fabs(slo) < 1e-8 )
        {
            if( fabs(sig) > 1e-4 )
                search_from = Infinity;
        }
        else
            search_from -= sig / slo;
    } while( old_search_from < search_from && search_from < next_incontinuity);

    if( debug_domain )
        cout << "leave domain_floor search_from=" << search_from << endl;
    return search_from;
}

real domain_ceil_time( vector<valdomain> &sig_domain, ntime_t next_incontinuity )
{
    if( debug_domain )
        cout << "enter domain_ceil" << endl;
    ntime_t search_from = 0;
    ntime_t old_search_from;

    do {
        if( debug_domain )
            cout << "domain_ceil search_from=" << search_from << endl;
        real sig = 0.0, slo = 0.0;
        for( vector<valdomain>::iterator i = sig_domain.begin(); i != sig_domain.end(); i++ )
        {
            if( debug_domain )
                cout << "domain " << *i << endl;
            real nsig = i->value + i->upslope * search_from;
            if( nsig >= i->ceil ) {
                sig += i->ceil;
            }
            else {
                sig += nsig; 
                slo += i->upslope;
            }
        }
        old_search_from = search_from;
        if( debug_domain )
            cout << "domain_ceil sig=" << sig << ", slo=" << slo << endl;
        if( fabs(slo) < 1e-8 )
        {
            if( fabs(sig) > 1e-4 )
                search_from = Infinity;
        }
        else
            search_from -= sig / slo;
    } while( old_search_from < search_from && search_from < next_incontinuity);

    if( debug_domain )
        cout << "leave domain_ceil search_from=" << search_from << endl;
    return search_from;
}

template <bool debug>
void tneuron_ext<debug>::pulseArrive(tscheduler &scheduler, ntime_t current_time, real pulse_level)
{
    real curval = pulses->getValue(current_time);

    pulses->setParam(curval + pulse_level, current_time);

    if( getDeb() )
    cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
        << " Pulse arrived, signal level " << pulse_level << ", new pulse " << pulses->getValue(current_time) << endl; 

    scheduleFire(scheduler, current_time);
}

template <bool debug>
void tneuron_ext<debug>::pulseArrive(tscheduler &scheduler, ntime_t current_time, message_base *base)
{
//  real curval = pulses->getValue(current_time);

    if( getDeb() )
    cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
        << " Pulse arrived, message=" << base->getMessageId() << flush;
    bool ret = sendMessage(current_time, base);
    if( getDeb() )
        cout << " result=" << ret << endl;

    scheduleFire(scheduler, current_time);
}

template <bool debug>
void tneuron_ext<debug>::pulseArrive(tscheduler &scheduler, ntime_t current_time, func_base *base)
{
//  real curval = pulses->getValue(current_time);

    if( getDeb() )
    cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
        << " Pulse arrived, func=" << base->getDescription() << flush;
    addExt(base);

    scheduleFire(scheduler, current_time);
}

template <bool debug>
void tneuron_ext<debug>::scheduleFire(tscheduler &scheduler, ntime_t current_time, bool /*resched*/)
{
    prof<1> profobj("scheduleFire");
//  if( getDeb() ) cout << "enter tneuron_ext::schedulefire" << endl;
    
    real signal = 0.0;
    real deriv1st = 0.0;
    real deriv2nd = 0.0;
    real next_incontinuity = Infinity;
    ntime_t next_known_pulse;

    int extssize = exts.size();
    vector<valdomain> sig_domain(extssize);
    vector<valdomain> d1st_domain(extssize);
    vector<valdomain> d2nd_domain(extssize);

    last_simulate = current_time;
    real maxgrad = 0.0;

    if( use_next_known_pulse ) { 
        prof<2> profobj("tneuron_ext::0 calcval");
        while( !queue()->empty() && &(queue()->top().getAction()) == this )
        {
            totalrescheduled++;
            queue()->pop();
        }
        if( queue()->empty() )
            next_known_pulse = Infinity;
        else
            next_known_pulse = queue()->top().getTime();

        vector<valdomain>::iterator sdi = sig_domain.begin();
        for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
        {
            real v = (*i)->getValue(current_time);
            sdi++->value = v;
            signal += v;

//          cout << "v=" << v << ", m=" << (*i)->getMaxGradient(current_time) << ", desc= "<< (*i)->getDescription() << endl;

            maxgrad += (*i)->getMaxGradient(current_time);
        }
    } 
    else {
        next_known_pulse = Infinity;
        vector<valdomain>::iterator sdi = sig_domain.begin();
        for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
        {
            real v = (*i)->getValue(current_time);
            valdomain vd; vd.value = v;
            sdi++->value = v;
            signal += v;
    //cout << "v=" << v << ", desc= "<< (*i)->getDescription() << endl;
        }
    }

    if( signal >= - 1e-8 )
    {
        if( getDeb() )
            cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
                << " Fire, signal level " << signal << endl; 
        fire(scheduler, current_time);
    }
    else
    {
        if( use_next_known_pulse ) { 
            if( signal < - maxgrad * (next_known_pulse - current_time) )
            {
                if( getDeb() )
                    cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
                        << " Stay, signal level " << signal << ", filtered by max gradient " << maxgrad << ", next " << next_known_pulse << endl;
                totalfiltered_maxgrad++;
                return;
            }
        }
        {
        prof<2> profobj("tneuron_ext::1 sigdomain");
        vector<valdomain>::iterator sdi = sig_domain.begin();
        for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
        {
    //      if( getDeb() )
    //          cout << "func_base " << (*i)->getDescription() << endl;
            while( (*i)->shouldDelete(current_time) )
            {
                prof<2> profobj("delete");
                delete *i;
                exts.erase(i);
                sig_domain.erase(sdi);
                if( i == exts.end() )
                    goto EXIT_DOM_LOOP;
            }

    //      if( sdi == sig_domain.end() ) cout << "Why? " << sig_domain.size() << " / " << exts.size() << endl;
            (*i)->getValueDomain(current_time, sdi->upslope, sdi->ceil, sdi->downslope, sdi->floor);
            sdi++;
        }
        EXIT_DOM_LOOP: ;
        }
        real value_bound;
        {
        prof<2> profobj("tneuron_ext::2 filter");
        value_bound = domain_ceil_time(sig_domain, next_incontinuity - current_time);

        if( value_bound + current_time >= next_incontinuity && next_incontinuity < Infinity )
        {
            if( getDeb() )
                cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
                    << " Stay, signal level " << signal << ", Reschedule at next incontinuity " << next_incontinuity << endl; 
            setLoopBack(scheduler, next_incontinuity);
            totalfiltered_incontinuity++;
            return;
        }
        if( use_next_known_pulse && value_bound + current_time >= next_known_pulse && next_known_pulse < Infinity )
        {
            if( getDeb() )
                cout << setw(9) << setiosflags(ios::fixed) << setprecision(debug_precision) << current_time << " " << name 
                    << " Stay, signal level " << signal << ", wait for the next known pulse at " << next_known_pulse << endl; 
            totalfiltered_nextpulse++;
            return;
        }
        }
        real next_bound, bound_val;
        {
        prof<2> profobj("tneuron_ext::3 bounder");
        if( iporder >= 1 )
        {
            vector<valdomain>::iterator d1d = d1st_domain.begin();
            for( vector<func_base *>::iterator i = exts.begin(); i != exts.end(); i++ )
            {
                real v = (*i)->get1stDeriv(current_time);
                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

---