package cnslab.cnsnetwork;
import java.util.*;
import cnslab.cnsmath.*;
import edu.jhu.mb.ernst.model.Synapse;
import edu.jhu.mb.ernst.util.slot.Slot;
/***********************************************************************
* MN model with STDP implemented.
*
* See Mihalas and Niebur 2009 paper.
*
* @version
* $Date: 2012-08-04 20:43:22 +0200 (Sat, 04 Aug 2012) $
* $Rev: 104 $
* $Author: croft $
* @author
* Yi Dong
* @author
* David Wallace Croft
* @author
* Jeremy Cohen
***********************************************************************/
public final class VSICLIFNeuronV2
implements Neuron
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
{
// public Axon axon;
// private double[] curr;//0 is exc current, 1 is inhib current
// private double memV;
private double timeOfNextFire;
private double timeOfLastUpdate;
/** whether the neuron is recordable or not */
private boolean record;
/** table which stores presynaptic event time and weight */
private Map<Synapse, TimeWeight> histTable;
//private double Alpha_LTD=1.0;//weight of strenght of LTD
//private double Alpha_LTP=1.0;//weight of strenght of LTP
//private double K_LTD=1/0.020; //LTD inverse time constant
//private double K_LTP=1/0.005; //LTP inverse time constant
//parameters were moved to VSICLIFNouronParaV2
// public static PrintStream p;
// public static boolean print=false;
/** minimum rising time for neuron to fire */
@Deprecated
public static double MINRISETIME = 1e-4;
/** time accuracy */
public static double tEPS = 1e-12;
/** the criteria for removing synapses, if a synaptic channel does not
get an input for a long time such that its current decays below rEPS,
it is removed from computations of the voltage and future updates up
until it gets an input */
public static double rEPS = 1e-22;
/** initialization of the maxWeight which was used in a previous method
in doesfire */
public static double maxWeight = 0;
/** the cost associated with neuron update,see D'hane's paper */
public static double cost_Update = 10;
/** the cost associated with queue schedule. see D'hane's paper */
public static double cost_Schedule = 1;
/** total cost, the costs associated to putting an temporary spike time
in */
public static double cost;
/** average time interval a neuron received spikes */
public double tAvg;
/** last time the neuron received a spike */
public double lastInputTime;
/** last time the neuron fired, negative means no last spike */
public double lastFireTime = -1.0;
/** clamp the membrane voltage during absolute refractory period */
public double clamp_v;
/** clamp the threshold during absolute refractory period */
public double clamp_t;
/** judge whether the neuron fires or not */
public boolean fire;
/** long imposes a constraint on the maximum number of hosts to be 64 */
public long tHost;
/** Neuron parameters */
public VSICLIFNeuronParaV2 para;
/** linked list for state variables */
public TreeSet<State> state;
/** pointer to the state variables; the first two states correspond to
the g and b terms respectively, followed by the spike-induced current
terms */
public State [ ] sta_p;
////////////////////////////////////////////////////////////////////////
// inner classes
////////////////////////////////////////////////////////////////////////
public final class State
implements Comparable<State>
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
{
public double time;
public double value;
public byte id;
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
public State (
final double time,
final double value,
final int id )
////////////////////////////////////////////////////////////////////////
{
this.time = time;
this.value = value;
this.id = ( byte ) id;
}
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
@Override
public int compareTo ( final State arg0 )
////////////////////////////////////////////////////////////////////////
{
if ( para.allDecays [ this.id ] < para.allDecays [ arg0.id ] )
{
return 1;
}
else if ( para.allDecays [ this.id ] > para.allDecays [ arg0.id ] )
{
return -1;
}
else
{
if ( this.id < arg0.id )
{
return -1;
}
else if ( this.id > arg0.id )
{
return 1;
}
else
{
return 0;
}
}
}
@Override
public String toString ( )
////////////////////////////////////////////////////////////////////////
{
return
"time:" + time
+ " value:" + value
+ " id:" + id
+ " decay:" + para.allDecays [ id ];
}
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
}
////////////////////////////////////////////////////////////////////////
// constructor methods
////////////////////////////////////////////////////////////////////////
public VSICLIFNeuronV2 ( final VSICLIFNeuronParaV2 para )
////////////////////////////////////////////////////////////////////////
{
this.para = para;
if ( maxWeight
< para.CAP * ( para.THRESHOLD - para.VRESET ) / para.MINRISETIME )
{
// not used
maxWeight
= para.CAP * ( para.THRESHOLD - para.VRESET ) / para.MINRISETIME;
}
}
////////////////////////////////////////////////////////////////////////
// interface Neuron accessor methods
////////////////////////////////////////////////////////////////////////
/*
public Axon getAxon()
{
return this.axon;
}
*/
@Override
public double [ ] getCurr ( final double currTime )
////////////////////////////////////////////////////////////////////////
{
final double [ ] out = new double [ para.DECAY.length ];
for ( int a = 0; a < para.DECAY.length; a++ )
{
if ( sta_p [ a + 2 ] == null )
{
out [ a ] = 0;
}
else
{
out [ a ] = sta_p [ a + 2 ].value * Math.exp (
-( currTime - sta_p [ a + 2 ].time ) * para.DECAY [ a ] )
* para.GLCAPDECAY [ a ] / ( 1 - para.ABDECAY [ a ] );
}
}
return out;
}
@Override
public double getMemV ( final double currTime )
////////////////////////////////////////////////////////////////////////
{
if ( currTime < timeOfLastUpdate )
{
// if input comes inside the refractory period, memV not changed
return memVoltage ( timeOfLastUpdate );
}
else
{
return memVoltage (currTime);
}
}
@Override
public boolean getRecord ( )
////////////////////////////////////////////////////////////////////////
{
return this.record;
}
@Override
public long getTHost ( )
////////////////////////////////////////////////////////////////////////
{
return tHost;
}
@Override
public double getTimeOfNextFire ( )
////////////////////////////////////////////////////////////////////////
{
return this.timeOfNextFire;
}
@Override
public boolean isSensory ( )
////////////////////////////////////////////////////////////////////////
{
return false;
}
@Override
public boolean realFire ( )
////////////////////////////////////////////////////////////////////////
{
//checks if the inserted spike is true,
// or if it a temporary placeholder
return fire;
}
////////////////////////////////////////////////////////////////////////
// interface Neuron mutator methods
////////////////////////////////////////////////////////////////////////
@Override
public void setRecord ( final boolean record )
////////////////////////////////////////////////////////////////////////
{
this.record = record;
}
@Override
public void setTHost ( final long id )
////////////////////////////////////////////////////////////////////////
{
this.tHost = id;
}
@Override
public void setTimeOfNextFire ( final double timeOfNextFire )
////////////////////////////////////////////////////////////////////////
{
this.timeOfNextFire = timeOfNextFire;
}
////////////////////////////////////////////////////////////////////////
// interface Neuron lifecycle methods
////////////////////////////////////////////////////////////////////////
@Override
public void init (
final int expid,
final int trialid,
final Seed idum,
final Network net,
final int id )
////////////////////////////////////////////////////////////////////////
{
// initialization all the initial parameters.
cost = Math.log ( cost_Schedule / cost_Update + 1.0 );
this.tAvg = 0.005; // initial value is 5 ms;
this.lastInputTime = 0.0;
this.timeOfLastUpdate = 0.0;
this.timeOfNextFire = -1;
double memV = para.ini_mem + para.ini_memVar * Cnsran.ran2 ( idum );
double instantaneousThreshold = para.ini_threshold;
double [ ] curr = new double [ para.ini_curr.length ];
for ( int a = 0; a < curr.length; a++ )
curr [ a ] = para.ini_curr [ a ];
para.CON_v = para.IEXT / para.GL;
para.CON = para.IEXT / para.GL - ( para.THRESHOLD - para.VREST
+ para.A / para.B * para.IEXT / para.GL );
// table history presynaptic event time and relative weight
histTable = new HashMap<Synapse,TimeWeight> ( );
// last time the neuron fired
this.lastFireTime = -1.0;
//state variables setup
//create linked table
state = new TreeSet<State> ( );
sta_p = new State [ para.DECAY.length + 2 ];
// First state is Ng term;
sta_p [ 0 ] = new State ( 0, getNG ( memV, curr ), 0 );
state.add ( sta_p [ 0 ] );
sta_p [ 1 ] = new State (
0,
getNB ( instantaneousThreshold, memV, curr ),
1 );
state.add ( sta_p [ 1 ] );
for ( int i = 2; i < para.DECAY.length + 2; i++ )
{
sta_p [ i ] = new State ( 0, getCurr ( i - 2, curr ), i );
state.add ( sta_p [ i ] );
}
/*
//System.out.println("starting");
Iterator<State> iter= state.iterator();
while(iter.hasNext())
{
//System.out.println(iter.next());
}
//System.out.println("ending");
*/
double dt, vc;
double de = safeD ( 0.0 );
double nextTime = ( -( memV - instantaneousThreshold ) / de );
dt = nextTime;
// if(nowV-nowT>0)return -1.0;
while ( de > 0 && nextTime < 1.0
&& !( dt > 0 ? dt < tEPS : -dt < tEPS ) )
{
de = safeD ( nextTime );
vc = memStheta ( nextTime );
dt = ( -vc / de );
// //System.out.print("n:"+nextTime+" s:"+(-vc)+"/"+de+"="
// // +dt+" Etime:"+(time+nextTime));
nextTime += dt;
}
if ( de <0 || nextTime > 1.0 )
{
fire = false;
nextTime = -1.0;
}
else if ( ( dt > 0 ? dt < tEPS : -dt < tEPS ) )
{
fire = true;
}
else
{
//System.out.println("error condition in init");
}
if ( nextTime < 0 )
{
// neuron will not fire at all
}
else
{
final Slot<FireEvent> fireEventSlot = net.getFireEventSlot ( );
// nextTime = (timeOfLastUpdate>0.0)
// ?(timeOfLastUpdate-0.0-Math.log(nextTime)/para.K)
// :(-1.0*Math.log(nextTime)/para.K);
if ( net.getClass ( ).getName().equals (
"cnslab.cnsnetwork.ANetwork" ) )
{
fireEventSlot.offer (
new AFireEvent (
id,
nextTime,
net.info.idIndex,
( ( ANetwork ) net ).aData.getId ( ( ANetwork ) net ) ) );
}
else if ( net.getClass ( ).getName ( ).equals (
"cnslab.cnsnetwork.Network" ) )
{
fireEventSlot.offer ( new FireEvent ( id, nextTime ) );
}
else
{
throw new RuntimeException (
"Other Network Class doesn't exist" );
}
timeOfNextFire = nextTime;
}
}
@Override
public double updateFire ( )
////////////////////////////////////////////////////////////////////////
{
// System.out.println("fire at:"+timeOfNextFire);
//update state variables containing membrane voltage and threshold
//because of the resetting.
double baseline;
double nowV, nowT;
if ( fire ) // for real fire
{
//Iterator<Synapse,WeightTime> it = histTable.iterator();
for ( final Map.Entry<Synapse,TimeWeight> entry
: histTable.entrySet ( ) )
{
final Synapse syn = entry.getKey ( );
// get the relative weight
final TimeWeight tw = entry.getValue();
// channel 0 has STDP
if ( syn.getType ( ) == 0 )
{
// LTP only for close spikes
if ( lastFireTime < tw.time )
{
//update the weight
tw.weight = tw.weight
* ( 1 + para.Alpha_LTP
* Math.exp ( -para.K_LTP
* ( timeOfNextFire-tw.time ) ) );
}
}
}
// Store the old neuron fire time.
lastFireTime = timeOfNextFire;
// backup current threshold
final double backCT
= memVoltage ( timeOfNextFire ) - memStheta ( timeOfNextFire );
// G term update
sta_p [ 0 ].value *= Math.exp (
-( timeOfNextFire - sta_p [ 0 ].time + para.ABSREF )
* para.allDecays [ 0 ] );
sta_p [ 0 ].time = timeOfNextFire + para.ABSREF;
// B term update
sta_p [ 1 ].value *= Math.exp (
-( timeOfNextFire-sta_p [ 1 ].time + para.ABSREF )
* para.allDecays [ 1 ] );
sta_p [ 1 ].time = timeOfNextFire + para.ABSREF;
//membrane voltage reset
//future memVoltage
final double
cVoltage = memVoltage ( timeOfNextFire + para.ABSREF );
// change memStheta(timeOfNextFire) to 0 later on may be beneficial
// double cThreshold = cVoltage - memStheta(timeOfNextFire);
// change memStheta(timeOfNextFire) to 0 later on may be beneficial
// double cThreshold
// = memVoltage(timeOfNextFire) - memStheta(timeOfNextFire);
// System.out.println (
// "abT:"+backCT+" memV:"+memVoltage(timeOfNextFire));
// System.out.println ( "abT2:"
// +(memVoltage(timeOfNextFire) - memStheta(timeOfNextFire)));
// membrane voltage offset
final double vOffset = para.VRESET - cVoltage;
// threshold in future
final double
cdiff = ( cVoltage - memStheta ( timeOfNextFire + para.ABSREF ) );
// threshold offset
final double tOffset
= Math.max (
para.RRESET,
backCT + para.THRESHOLDADD )
- cdiff;
sta_p [ 0 ].value += ( 1 - para.ABGLCAP ) * vOffset;
sta_p [ 1 ].value += para.ABGLCAP * vOffset - tOffset;
// System.out.println("after fire voltage:"
// +memVoltage(timeOfNextFire+para.ABSREF)
// +" v-t:"+memStheta(timeOfNextFire+para.ABSREF));
// if(clamp_v==Double.MAX_VALUE) //if clamp_v is not set;
clamp_v = para.VRESET;
clamp_t = Math.max ( para.RRESET, backCT + para.THRESHOLDADD );
// this section is for spike triggered current
for ( int a = 0; a < para.DECAY.length; a++ )
{
// update the currents if spike trigger current is non zero
if ( para.IRATIO [ a ] != 1.0 || para.IADD [ a ] != 0 )
{
// System.out.println("added");
// update current to fire time;
if ( sta_p [ a + 2 ] == null ) //activate it
{
//System.out.println("act");
sta_p [ a + 2 ] = new State (
timeOfNextFire + para.ABSREF,
0,
a + 2 );
state.add ( sta_p [ a + 2 ] );
}
else
{
sta_p [ a + 2 ].value *= Math.exp (
-( timeOfNextFire - sta_p [ a + 2 ].time + para.ABSREF )
* para.allDecays [ a + 2 ] );
sta_p [ a + 2 ].time = timeOfNextFire + para.ABSREF;
}
// current at timeofnextfire
final double currNow = sta_p [ a + 2 ].value
* para.GLCAPDECAY [ a ]
/ ( 1 - para.ABDECAY [ a ] );
final double currAfter = ( para.IRATIO [ a ] * currNow
* Math.exp ( para.ABSREF * para.allDecays [ a + 2 ] )
+ para.IADD [ a ] )
* Math.exp ( -para.ABSREF * para.allDecays [ a + 2 ] );
final double changeCurr = currAfter - currNow;
// System.out.print("after:"+ currAfter+" now:"+currNow+" ");
sta_p[0].value
+= changeCurr / para.GLCAPDECAY [ a ] * ( para.ABGLCAP - 1 );
sta_p[1].value
-= changeCurr / para.BCAPDECAY [ a ] * para.ABGLCAP;
sta_p[a+2].value
+= ( changeCurr / para.GLCAPDECAY [ a ] )
* ( 1 - para.ABDECAY [ a ] );
/*
sta_p[0].value
+= input.weight/para.GLCAPDECAY[a]*(para.ABGLCAP-1);
sta_p[1].value
-= input.weight/para.BCAPDECAY[a]*para.ABGLCAP;
sta_p[a+2].value
+= (changeCurr/para.GLCAPDECAY[a])*(1 - para.ABDECAY[a]);
*/
// System.out.print ( " bycalcu:" + sta_p [ a + 2 ].value
// * para.GLCAPDECAY [ a ] / ( 1 - para.ABDECAY [ a ] )
// + " " );
}
}
// System.out.println ( "after fire curr:"
// + getCurr ( timeOfNextFire + para.ABSREF ) [ 0 ]
// +" "+getCurr(timeOfNextFire+para.ABSREF)[1]);
timeOfLastUpdate = timeOfNextFire+para.ABSREF;
baseline = timeOfNextFire+para.ABSREF;
nowV = para.VRESET;
nowT = Math.max ( para.RRESET, backCT + para.THRESHOLDADD );
}
else //neuron is partial updated.
{
baseline = timeOfNextFire;
nowV = memStheta(baseline);
nowT = 0;
}
double dt, vc;
double de = safeD ( baseline );
double nextTime = ( -( nowV - nowT ) ) / de;
dt = nextTime;
while ( de > 0 && nextTime < cost * tAvg )
{
de = safeD ( baseline + nextTime );
vc = memStheta ( baseline + nextTime );
dt = ( -vc / de );
if ( dt > 0 ? dt < tEPS : -dt < tEPS )
{
fire = true;
//no spike within absolute refractory period.
//if(nextTime < para.ABSREF) return -1.0;
//System.out.println(
// "Again:"+(nextTime+baseline - timeOfNextFire));
return nextTime + baseline - timeOfNextFire;
}
nextTime += dt;
}
if ( de < 0 || nextTime > 1.0 )
{
fire = false;
return -1.0;
}
else
{
fire = false;
return nextTime + baseline - timeOfNextFire;
}
}
@Override
public double updateInput (
final double time,
final Synapse input )
////////////////////////////////////////////////////////////////////////
{
if ( histTable.containsKey ( input ) )
{
// if this synapse is already in the table, just update its value
TimeWeight tw = histTable.get ( input );
if ( input.getType ( ) == 0 ) //channel 0 has STDP
{
if ( lastFireTime > tw.time ) // LTD only for close spikes
{
// update weight if neuron fired before
tw.weight = tw.weight * ( 1 - para.Alpha_LTD * Math.exp (
-para.K_LTD * ( time - lastFireTime ) ) );
// System.out.println("LTD:"+tw.weight);
}
}
tw.time = time; // update the time
}
else // if this synapse did not fire before, add it to the table
{
TimeWeight tw = new TimeWeight ( time, 1.0 );
if ( input.getType ( ) == 0 ) //channel 0 has STDP
{
if ( lastFireTime > 0.0 )
{
// update weight if neuron fired before
tw.weight = tw.weight * ( 1 - para.Alpha_LTD * Math.exp (
-para.K_LTD * ( time - lastFireTime ) ) );
}
}
// put the default weight into the history table
histTable.put ( input, tw );
}
// System.out.println ( "get input at:" + time + " syn:" + input
// + " membrane voltage is:" + memVoltage ( time ) );
// System.out.println ( "beg: t:" + time + " s:" + input
// + " mV:" + memVoltage ( time ) + " mtV:" + memStheta ( time ) );
//update the tAvg; look for 5 points in the past
tAvg = tAvg * 0.8 + ( time - lastInputTime ) *0.2;
lastInputTime = time;
// System.out.println ( "before updateInput: memV: " + memV
// + " instantThresh: " + instantaneousThreshold + " curr0: "
// + curr[0] + " curr1: " + curr[1]);
//
// p.println("updateInput"+time);
// p.flush();
/*
double [] oldcurr = new double [curr.length];
oldcurr[0] = curr[0];
oldcurr[1] = curr[1];
double oldMem=memV ;
double VACC = 1E-9; // accuracy in voltage for spike to be generated
*/
// System.out.println ( "0t:" + time + " s:" + input
// + " mV:" + memVoltage ( time ) + " mtV:" + memStheta ( time )
// + " th:" + ( memVoltage ( time ) - memStheta ( time ) ) );
// update state 0,1 to current time;
sta_p [ 0 ].value *= Math.exp (
-( time - sta_p [ 0 ].time ) * para.allDecays [ 0 ] );
// System.out.println ( "ng decay:" + Math.exp (
// -( time - sta_p [ 0 ].time ) * para.allDecays [ 0 ] ) );
sta_p [ 0 ].time = time;
sta_p [ 1 ].value *= Math.exp (
-( time - sta_p [ 1 ].time ) * para.allDecays [ 1 ] );
sta_p [ 1 ].time = time;
// System.out.println ( "mid1: t:" + time + " s:" + input + " mV:"
// + memVoltage ( time ) + " mtV:" + memStheta ( time ) );
//check whether the input is deleted
if ( sta_p [ input.getType ( ) + 2 ] == null ) //activate it
{
//System.out.println("act");
sta_p [ input.getType ( ) + 2 ] = new State ( time, 0, input.getType ( ) + 2 );
state.add ( sta_p [ input.getType ( ) + 2 ] );
}
else
{
sta_p [ input.getType ( ) + 2 ].value *= Math.exp ( -( time - sta_p [
input.getType ( ) + 2].time ) * para.allDecays [ input.getType ( ) + 2 ] );
sta_p[input.getType ( )+2].time = time;
}
// System.out.println ( "mid2: t:" + time + " s:" + input + " mV:"
// + memVoltage ( time ) + " mtV:" + memStheta ( time ) );
// add the extra state values (0,2,this synapse) by getting input
// spike.
sta_p [ 0 ].value += input.getWeight ( ) * ( histTable.get ( input ).weight )
/ para.GLCAPDECAY [ input.getType ( ) ] * ( para.ABGLCAP - 1 );
sta_p [ 1 ].value -= input.getWeight ( ) * ( histTable.get ( input ).weight )
/ para.BCAPDECAY [ input.getType ( ) ] * para.ABGLCAP;
sta_p [ input.getType ( ) + 2 ].value
+= ( input.getWeight ( ) * ( histTable.get ( input ).weight )
/ para.GLCAPDECAY [ input.getType ( ) ] )
* ( 1 - para.ABDECAY [ input.getType ( ) ] );
double nextTime;
double baseline;
double nowV = memVoltage(time);
double nowT = nowV - memStheta(time);
if ( time < timeOfLastUpdate )
{
// set the offset of the membraneVoltage to make sure at the
// timeOfLastUpdate the membraneVoltage is not changed at all
// CON_v = CON_v - (memStheta(timeOfLastUpdate) - nowMem);
/*
if(clamp_v==Double.MAX_VALUE) //if clamp_v is not set;
{
clamp_v= nowV;
clamp_t= nowT;
}
*/
double vFuture = memVoltage ( timeOfLastUpdate );
double voltageOffset = ( vFuture - clamp_v );
double tOffset = vFuture - memStheta ( timeOfLastUpdate ) - clamp_t;
sta_p [ 0 ].value += -( 1 - para.ABGLCAP ) * voltageOffset
* Math.exp ( ( timeOfLastUpdate - time ) * para.allDecays [ 0 ] );
sta_p [ 1 ].value += -( para.ABGLCAP * voltageOffset - tOffset )
* Math.exp ( ( timeOfLastUpdate - time ) * para.allDecays [ 1 ] );
// System.out.println ( "nowV:" + nowV
// + " futV:" + memVoltage ( timeOfLastUpdate )
// + " clamp_v:" + clamp_v );
// System.out.println ( "nowT:" + nowT
// + " futT:" + ( memVoltage ( timeOfLastUpdate )
// - memStheta ( timeOfLastUpdate ) )
// + " clamp_6:" + clamp_t );
// sta_p[0].value += (1-para.ABGLCAP)*vOffset;
// sta_p[1].value += para.ABGLCAP*vOffset - tOffset;
nowV = clamp_v;
nowT = clamp_t;
baseline = timeOfLastUpdate;
}
else
{
clamp_v=Double.MAX_VALUE;
clamp_t=Double.MAX_VALUE;
//if(input.to==1)
{
//System.out.println("t:"+time+ " s:"+input+ " mV:"
//+ memVoltage(time) +" mtV:"+memStheta(time)+ " th:"
//+(memVoltage(time)-memStheta(time))+ " last:"+timeOfLastUpdate);
// //System.out.println("t:"+time+ " s:"+input+ " mV:"+memV
// +" mtV:"+ (memV- instantaneousThreshold)+ " th:"
// +instantaneousThreshold);
// double [] out= getCurr(time);
//System.out.println("c0:"+out[0]+ " c1:"+out[1]);
}
timeOfLastUpdate=time;
baseline = time;
}
double dt, vc;
double de = safeD ( baseline );
nextTime = ( -( nowV - nowT ) / de );
dt = nextTime;
// if(nowV-nowT>0)return -1.0;
while ( de > 0 && nextTime < cost * tAvg )
{
de = safeD ( baseline + nextTime );
vc = memStheta ( baseline + nextTime );
dt = ( -vc / de );
// //System.out.print("n:"+nextTime+" s:"+(-vc)+"/"+de+"="+dt
// +" Etime:"+(time+nextTime));
if ( dt > 0 ? dt < tEPS : -dt < tEPS )
{
fire = true;
// //System.out.println();
// if(nextTime < timeOfLastUpdate - time)
// {
// //System.out.println("fail to fire nextTime:"+nextTime
// //+" smaller than:"+(timeOfLastUpdate - time));
// return -1.0; //no spike within absolute refractory period.
// }
//System.out.println("N :"+input.to+" fire:"
//+(time>timeOfLastUpdate?nextTime
//:nextTime+timeOfLastUpdate-time)
//+ " memV:"+memStheta(baseline+nextTime));
return ( time > timeOfLastUpdate
? nextTime : nextTime + timeOfLastUpdate - time );
}
nextTime += dt;
}
if( de < 0 || nextTime > 1.0 )
{
// input is not strong enough to make neuron fire
fire = false;
return -1.0;
}
else
{
// neuron dot fire yet, but eventually will fire in future,
// do a partial update.
fire = false;
return ( time > timeOfLastUpdate
? nextTime : nextTime + timeOfLastUpdate - time );
}
// //System.out.println();
// //System.out.println("fail to fire dt:"+dt+ " nextTime:"
// // +nextTime + " test:"+memStheta(time+0.10316241162745285)
// // + " deri:"+de+ " safeD:"+safeD(time+0.10316241162745285));
//return -1.0;
}
////////////////////////////////////////////////////////////////////////
// accessor methods
////////////////////////////////////////////////////////////////////////
/***********************************************************************
* get the current which is the last updated current
***********************************************************************/
public double [ ] getCurr ( )
////////////////////////////////////////////////////////////////////////
{
final double [ ] out = new double [ para.DECAY.length ];
for ( int a = 0; a < para.DECAY.length; a++ )
{
if ( sta_p [ a + 2 ] == null )
{
out [ a ] = 0;
}
else
{
out [ a ] = sta_p [ a + 2 ].value * para.GLCAPDECAY [ a ]
/ ( 1 - para.ABDECAY [ a ] );
}
}
return out;
}
/***********************************************************************
* Nj terms:
***********************************************************************/
public double getCurr (
final int a,
double [ ] curr )
////////////////////////////////////////////////////////////////////////
{
//Nj terms:
return ( curr [ a ] / para.GLCAPDECAY [ a ] )
* ( 1 - para.ABDECAY [ a ] );
}
/***********************************************************************
* return nb term;
***********************************************************************/
public double getNB (
final double instantaneousThreshold,
final double memV,
final double [ ] curr ) //return nb term;
////////////////////////////////////////////////////////////////////////
{
double nbTerm
= para.VREST - memV + para.IEXT / ( para.B * para.CAP );
for ( int a = 0; a < curr.length; a++ )
nbTerm += curr [ a ] / para.BCAPDECAY [ a ];
// replace b*cap for para.GL
nbTerm *= para.ABGLCAP;
nbTerm += instantaneousThreshold - para.THRESHOLD;
return -nbTerm;
}
/***********************************************************************
* Get ng term ( membrane voltage - ABGLCAP*membranevoltageNg )
***********************************************************************/
public double getNG (
final double memV,
final double [ ] curr )
////////////////////////////////////////////////////////////////////////
{
// return ng term ( membrane voltage - ABGLCAP*membranevoltageNg );
double ngTerm1 = memV - para.VREST - para.IEXT / para.GL;
for ( int a = 0; a < curr.length; a++ )
ngTerm1 -= curr [ a ] / para.GLCAPDECAY [ a ];
return ngTerm1 - ngTerm1 * para.ABGLCAP;
}
public double getSensoryWeight ( )
////////////////////////////////////////////////////////////////////////
{
throw new RuntimeException ( "This neuron type doesn't use the "
+ "Sensory Weight Functions!" );
}
public double getTimeOfLastUpdate ( )
////////////////////////////////////////////////////////////////////////
{
return this.timeOfLastUpdate;
}
/***********************************************************************
* get the membrane voltage
*
* @param t
* absolute time t
***********************************************************************/
public double memVoltage ( final double t )
////////////////////////////////////////////////////////////////////////
{
/*
double V=0.0;
Iterator<State> iter = state.iterator();
while(iter.hasNext())
{
State tmpState = iter.next();
if((int)tmpState.id != 1 ) // nb term will be ignored.
{
if(t-tmpState.time!=0)
{
if((int)tmpState.id == 0 ) // ng term
{
V += Math.exp(-(t-tmpState.time)
*para.allDecays[(int)tmpState.id])
**tmpState.value/(1-para.ABGLCAP);
// //System.out.print("ng:"+(tmpState.value/(1-para.ABGLCAP))
// //+ " de:"+Math.exp(-(t-tmpState.time)
// //*para.allDecays[(int)tmpState.id]));
}
else //current terms
{
V += Math.exp(-(t-tmpState.time)*para.allDecays[
(int)tmpState.id])*tmpState.value/(1 - para.ABDECAY[
(int)tmpState.id-2]);
// //System.out.print(" c:"+(tmpState.id-2)+" "+(
// //tmpState.value/(1 - para.ABDECAY[(int)tmpState.id-2]))
// //+ " de:"+Math.exp(-(t-tmpState.time)*para.allDecays[
// //(int)tmpState.id]));
}
}
else
{
if((int)tmpState.id == 0 ) // ng term
{
V += tmpState.value/(1-para.ABGLCAP);
// //System.out.print(" ng:"+(tmpState.value/(1-para.ABGLCAP))
// //+ " de:0");
}
else //current terms
{
V += tmpState.value/(1 - para.ABDECAY[(int)tmpState.id-2]);
// //System.out.print(" c:"+(tmpState.id-2)+" "
// //+(tmpState.value/(1 - para.ABDECAY[(int)tmpState.id-2]))
// //+" de:0");
}
}
}
}
// //System.out.print("V:"+V+" memV:"+(V+para.CON_v+para.VREST)+"\n");
// //System.out.println("V:"+V+" memV:"+(V+para.CON_v+para.VREST));
//System.out.println((V+para.CON_v+para.VREST));
*/
double V = 0;
double constant = 0;
Iterator<State> iter = state.iterator ( );
State first = iter.next ( );
State second;
double firstV, secondV;
if ( first.id == 1 ) //ignore nb term
{
if ( iter.hasNext ( ) )
{
first = iter.next ( );
}
else
{
return V + constant + para.CON_v + para.VREST;
}
}
if ( first.time == t ) // ignore no time change term
{
if ( first.id == 0 )
{
constant += first.value / ( 1 - para.ABGLCAP );
}
else
{
constant
+= first.value / ( 1 - para.ABDECAY [ ( int ) first.id - 2 ] );
}
if ( iter.hasNext ( ) )
{
first = iter.next ( );
}
else
{
return V + constant + para.CON_v + para.VREST;
}
}
if ( first.id == 1 ) //ignore nb term
{
if ( iter.hasNext ( ) )
{
first = iter.next ( );
}
else
{
return V + constant + para.CON_v + para.VREST;
}
}
if ( first.id == 0 )
{
firstV = first.value / ( 1 - para.ABGLCAP );
}
else
{
//System.out.println(first.id);
firstV
= first.value / ( 1 - para.ABDECAY [ ( int ) first.id - 2 ] );
}
V = firstV;
while ( iter.hasNext ( ) )
{
second = iter.next ( );
secondV = 0.0;
if ( second.time == t && second.id != 1 )
{
if ( second.id == 0 )
{
constant += second.value / ( 1 - para.ABGLCAP );
}
else
{
constant += second.value
/ ( 1 - para.ABDECAY [ ( int ) second.id - 2 ] );
}
if ( iter.hasNext ( ) )
{
second = iter.next ( );
}
else
{
break;
}
}
if ( second.id == 1 )
{
if ( iter.hasNext ( ) )
{
second = iter.next ( );
}
else
{
break;
}
}
if ( second.id == 0 )
{
secondV = second.value / ( 1 - para.ABGLCAP );
}
else
{
secondV = second.value
/ ( 1 - para.ABDECAY [ ( int ) second.id - 2 ] );
}
V = Math.exp ( -( ( para.allDecays [ ( int ) first.id ]
- para.allDecays [ ( int ) second.id ] ) * t
- ( para.allDecays [ ( int ) first.id ] * first.time
- para.allDecays [ ( int ) second.id ] * second.time ) ) ) * V
+ secondV;
first = second;
}
V = Math.exp ( -( para.allDecays [ ( int ) first.id ] * t
- para.allDecays [ ( int ) first.id ] * first.time ) ) * V;
//// System.out.println("V:"+(V+para.CON_v+para.VREST));
return V + para.CON_v + para.VREST + constant;
}
////////////////////////////////////////////////////////////////////////
// mutator methods
////////////////////////////////////////////////////////////////////////
/***********************************************************************
* set the membrane voltage
***********************************************************************/
public void setMemV ( double memV )
////////////////////////////////////////////////////////////////////////
{
return;
}
public void setTimeOfLastUpdate ( final double timeOfLastUpdate )
////////////////////////////////////////////////////////////////////////
{
this.timeOfLastUpdate = timeOfLastUpdate;
}
////////////////////////////////////////////////////////////////////////
// overridden Object methods
////////////////////////////////////////////////////////////////////////
@Override
public String toString ( )
////////////////////////////////////////////////////////////////////////
{
String tmp="";
// tmp = tmp + "MemVol:" + memV + "\n";
tmp=tmp+"Current:"+"\n";
// for(int i=0;i<curr.length;i++)
// {
// tmp=tmp+" curr"+i+" "+curr[i]+"\n";
// }
// tmp=tmp+axon;
return tmp;
}
////////////////////////////////////////////////////////////////////////
// miscellaneous methods
////////////////////////////////////////////////////////////////////////
/***********************************************************************
* Membrane voltage computate mem = S+\theta
*
* @param t
* absolute time t
***********************************************************************/
public double memStheta ( double t )
////////////////////////////////////////////////////////////////////////
{
/*
double V=0.0;
Iterator<State> iter = state.iterator();
while(iter.hasNext())
{
State tmpState = iter.next();
if(t-tmpState.time!=0)
{
V += Math.exp(-(t-tmpState.time)*para.allDecays[
(int)tmpState.id])*tmpState.value;
}
else
{
V += tmpState.value;
}
}
return V+para.CON;
*/
double V = 0, constant = 0.0;
Iterator<State> iter = state.iterator ( );
State first = iter.next ( );
State second;
double firstV, secondV;
if ( Math.abs ( first.value ) < rEPS && ( int ) first.id > 1 )
{
// delete the synapse if it is small
iter.remove ( );
sta_p [ ( int ) first.id ] = null;
}
// ignore no time change term
if ( first.time == t )
{
constant += first.value;
if ( iter.hasNext ( ) )
{
first = iter.next ( );
}
else
{
return V + constant + para.CON;
}
}
firstV = first.value;
V = firstV;
while ( iter.hasNext ( ) )
{
second = iter.next ( );
if ( Math.abs ( second.value ) < rEPS && ( int ) second.id > 1 )
{
// delete the synapse if it is small
iter.remove ( );
sta_p [ ( int ) second.id ] = null;
}
secondV = 0.0;
// ignore no time change term
if ( second.time == t )
{
constant += second.value;
if ( iter.hasNext ( ) )
{
second = iter.next ( );
}
else
{
break;
}
}
secondV = second.value;
V = Math.exp ( -( ( para.allDecays [ ( int ) first.id ]
- para.allDecays [ ( int ) second.id ] ) * t
- ( para.allDecays [ ( int ) first.id ] * first.time
- para.allDecays [ ( int ) second.id ] * second.time ) ) )
* V + secondV;
first = second;
}
V = Math.exp ( -( para.allDecays [ ( int ) first.id ] * t
- para.allDecays [ ( int ) first.id ] * first.time ) ) * V;
////System.out.println("V:"+(V+para.CON_v+para.VREST));
return V + para.CON + constant;
}
/***********************************************************************
* Safe derivative, good for Newton method, see D'hane paper
*
* @param t
* absolute time t
***********************************************************************/
public double safeD ( double t )
////////////////////////////////////////////////////////////////////////
{
// the smallest inverse decay of all negative state variables
double tauSafe = Double.MAX_VALUE;
//find the tauSafe
Iterator<State> iter = state.descendingIterator();
while(iter.hasNext())
{
State tmpState = iter.next();
// //System.out.println("value:"
// // +tmpState.value+" delay:"+para.allDecays[(int)tmpState.id]
// // + " id:"+(int)tmpState.id+ " time:"+tmpState.time);
if ( tmpState.value < 0 )
{
tauSafe = para.allDecays [ ( int ) tmpState.id ];
break;
}
}
////System.out.println("
// //System.out.println("safe:"+tauSafe);
//calculate the dV
//
double constant = 0;
double V = 0;
iter = state.iterator ( );
State first = iter.next ( );
State second;
double firstV, secondV;
if ( first.time == t )
{
// ignore no time change term
if ( first.value > 0 )
{
final double maxDecay = Math.min (
para.allDecays [ ( int ) first.id ],
tauSafe );
constant += -first.value * maxDecay;
}
else
{
constant += -first.value * para.allDecays [ ( int ) first.id ];
}
if ( iter.hasNext ( ) )
{
first = iter.next ( );
}
else
{
return V + constant;
}
}
if ( first.value > 0 )
{
final double maxDecay = Math.min (
para.allDecays [ ( int ) first.id ],
tauSafe );
firstV = -first.value * maxDecay;
}
else
{
firstV = -first.value * para.allDecays [ ( int ) first.id ];
}
V = firstV;
while ( iter.hasNext ( ) )
{
second = iter.next ( );
secondV = 0.0;
if( second.time == t)
{
if ( second.value > 0 )
{
final double maxDecay = Math.min (
para.allDecays [ ( int ) second.id ],
tauSafe );
constant += -second.value * maxDecay;
}
else
{
constant
+= -second.value
* para.allDecays [ ( int ) second.id ];
}
if ( iter.hasNext ( ) )
{
second = iter.next ( );
}
else
{
break;
}
}
if ( second.value > 0 )
{
final double maxDecay = Math.min (
para.allDecays [ ( int ) second.id ],
tauSafe );
secondV = -second.value * maxDecay;
}
else
{
secondV = -second.value * para.allDecays [ ( int ) second.id ];
}
V = Math.exp ( -( ( para.allDecays [ ( int ) first.id ]
- para.allDecays [ ( int ) second.id ] ) * t
- ( para.allDecays [ ( int ) first.id ] * first.time
- para.allDecays [ ( int ) second.id ] * second.time ) ) ) * V
+ secondV;
first = second;
}
V = Math.exp ( -( para.allDecays [ ( int ) first.id ] * t
- para.allDecays [ ( int ) first.id ] * first.time ) ) * V;
/*
double dV=0;
iter = state.iterator();
while(iter.hasNext())
{
State tmpState = iter.next();
if(tmpState.value>0)
{
double maxDecay = Math.min(
para.allDecays[(int)tmpState.id],tauSafe);
if(t-tmpState.time!=0)
{
// dV -= Math.exp(-(t-tmpState.time)*maxDecay)
// * tmpState.value*maxDecay;
dV -= Math.exp(-(t-tmpState.time)*para.allDecays[
(int)tmpState.id])*tmpState.value*maxDecay;
}
else
{
dV -= tmpState.value*maxDecay;
}
}
else
{
if(t-tmpState.time!=0)
{
dV -= Math.exp(-(t-tmpState.time)*para.allDecays[
(int)tmpState.id])*tmpState.value*para.allDecays[
(int)tmpState.id];
}
else
{
dV -= tmpState.value*para.allDecays[(int)tmpState.id];
}
}
}
*/
return V + constant;
}
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
}