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 13:43:22 -0500 (Sat, 04 Aug 2012) $
    *   $Rev: 104 $
    *   $Author: croft $
    * @author
    *   Yi Dong
    * @author
    *   David Wallace Croft, M.Sc.
    * @author
    *   Jeremy Cohen
    ***********************************************************************/
    public final class  MiNiNeuron
      implements Neuron
    ////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////
    {
    
    private double  timeOfNextFire;
    
    private double  timeOfLastUpdate;
    
    /** whether the neuron is recordable or not */
    private boolean  record;

    /** table for STDP that stores presynaptic event time and weight */
    private Map<Synapse, TimeWeight>  histTable;
    
    /** the epsilon of a final firing time calculation */
    public static double  tEPS = 1e-12;
    
    /** the criterion 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;

    /** the cost associated with one neuron update (see D’Haene, 2009) */
    public static double  cost_Update = 10;
    
    /** the cost associated with one queue schedule (see D'Haene, 2009) */
    public static double  cost_Schedule = 1;
    
    /** the expected cost of inserting a preliminary fire estimate into the 
     * event queue (as opposed to performing another NR iteration) */
    public static double  cost;

    /** running average of the time interval between input 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  clampVoltage;
    
    /** clamp the threshold during absolute refractory period */
    public double  clampThreshold;
    
    /** whether the next fire event is a final calculation,
     *  as opposed to a preliminary estimate (see D'Haene, 2009) */
    public boolean  fire; 

    /** long imposes a constraint on the maximum number of hosts to be 64 */
    public long  tHost;
    
    /** neuron parameters */
    public MiNiNeuronPara  para;

    /** linked list for state variables */
    public TreeSet<State>  state;
    
  /**
   * two-dimensional array pointing to the state variables. The first 
   * dimension corresponds to type -- NG, NB, NJ_SPIKE, and NJ_SYNAPSE, in 
   * that order. The second dimension corresponds to index. The NG and NB
   * arrays have only one element.
   */
    public State [ ] [ ] sta_p;
    

    ////////////////////////////////////////////////////////////////////////
    // inner classes
    ////////////////////////////////////////////////////////////////////////
    
    public final class  State
      implements Comparable<State>
    ////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////
    {
    /**
     * values of <code>type</code>  
     */
    public final static int 
                        NG         = 0,
                                NB         = 1,
                                NJ_SPIKE   = 2,
                                NJ_SYNAPSE = 3;
  
    /**
     * The time when this State was last updated.
     */
    public double  time;
   
    /**
     * The value of this State at the time when it was last updated.
     */
    public double  value;

    /**
     * The term in the thresholdDiff function that this State represents
     * (see Mihalas and Niebur, 2009, Equation 3.5.)
     */
    public int     type;
    
    /**
     * The index of this State among other States of the same type.
     */
    public int     index;

    ////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////

    public State (
      final double  time,
      final double  value,
      final int     type,
      final int     index)
    ////////////////////////////////////////////////////////////////////////
    {
      this.time  = time;

      this.value = value;

      this.type = type;
      
      this.index = index;
    }

    ////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////

    /**
     * Sorts in descending order of decay rate.
     */
    @Override
    public int  compareTo ( final State  arg0 )
    ////////////////////////////////////////////////////////////////////////
    {
        if ( para.allDecays [ this.type ] [ this.index ] <
                para.allDecays [ arg0.type ] [ arg0.index ])
        {
            return 1;
        }
        else if ( para.allDecays [ this.type ] [ this.index ] >
        para.allDecays [ arg0.type ] [ arg0.index ])
        {
            return -1;
        }

        if ( this.type < arg0.type )
        {
            return -1;
        }
        else if ( this.type > arg0.type )
        {
            return 1;
        }

        if (this.index < arg0.index)
        {
            return -1;
        }
        else if (this.index > arg0.index)
        {
            return 1;
        }

        else
        {
            return 0;
        }
    }

    @Override
    public String  toString ( )
    ////////////////////////////////////////////////////////////////////////
    {
      return
        "time:"   + time
        + " value:" + value
        + " type:"  + type
        + " index:" + index
        + " decay:" + para.allDecays [ type ] [ index ];
    }

    ////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////
    }

    ////////////////////////////////////////////////////////////////////////
    // constructor methods
    ////////////////////////////////////////////////////////////////////////
    
    public  MiNiNeuron ( final MiNiNeuronPara  para )
    ////////////////////////////////////////////////////////////////////////
    {
      this.para = para;
    }

    ////////////////////////////////////////////////////////////////////////
    // interface Neuron accessor methods
    ////////////////////////////////////////////////////////////////////////

    /**
     * Used mainly for intracellular recording.
     * 
     * @param currTime the current time.
     * @return the synaptic currents in this neuron at time currTime.
     */
    @Override
    public double [ ]  getCurr ( final double  currTime )
    ////////////////////////////////////////////////////////////////////////
    {
      final double [ ]  out = 
          new double [ sta_p [ State.NJ_SYNAPSE ].length ];
      
      for ( int  a = 0; a < sta_p [ State.NJ_SYNAPSE ].length; a++ )
      {
        State state = sta_p [ State.NJ_SYNAPSE ] [ a ];
        
        if ( state == null )
        {
          out [ a ] = 0;
        }
        else
        {
          // Convert NJ_SYNAPSE state to synapse current, and update to 
          // currTime.
          
          out [ a ] = state.value * Math.exp (
            -( currTime - state.time )
            * para.allDecays [ State.NJ_SYNAPSE ] [ a ] )
            / para.SYNAPSE_CURRENT_TO_SYNAPSE_NJ [ a ];
        }
      }
      
      return out;
    }    
    
    /**
     * Used mainly for intracellular recording.
     * 
     * @param currTime the current time.
     * @return the membrane voltage in this neuron at time currTime.
     */
    @Override
    public double  getMemV ( final double  currTime )
    ////////////////////////////////////////////////////////////////////////
    {
      if ( currTime < timeOfLastUpdate )
      {
        // If input comes inside the refractory period,
        // the membrane voltage remains unchanged.
        
        return membraneVoltage ( timeOfLastUpdate );
      }
      else
      {       
        return membraneVoltage (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 ( )
    ////////////////////////////////////////////////////////////////////////
    {
      // Whether the next fire event is a final calculation,
      // as opposed to a preliminary estimate (see D'Haene, 2009).
      
      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  initialThreshold = para.ini_threshold;
      
      // Pick the initial membrane voltage from a random, uniform spread.
      
      double  initialVoltage = 
          para.ini_mem + para.ini_memVar * Cnsran.ran2 ( idum );

      double [ ]  initialSpikeCurrents =
          new double [ para.ini_spike_curr.length ];
      
      for ( int  a = 0; a < initialSpikeCurrents.length; a++ )
        initialSpikeCurrents [ a ] = para.ini_spike_curr [ a ];
      
      // Set up the history table to moniter of synaptic activity for STDP
      
      histTable = new HashMap<Synapse,TimeWeight> ( );
      
      // Last time the neuron fired.
      
      this.lastFireTime = -1.0;
      
      // Set up state variables.

      // Create linked table for states.
      
      state = new TreeSet<State> ( );
      
      sta_p = new State [ 4 ] [ ];
      
      // Create the NG term.
      
      double ngTerm = para.NG_BASE;
      
      ngTerm += initialVoltage * para.VOLTAGE_TO_NG;
      
      for (int a = 0; a < initialSpikeCurrents.length; a ++)
      {
          ngTerm += 
              initialSpikeCurrents [ a ] * para.SPIKE_CURRENT_TO_NG [ a ];
      }      
      
      sta_p [ State.NG ] = new State [ ] { 
                  new State ( 0, ngTerm, State.NG, 0 ) 
              };
      
      state.add ( sta_p [ State.NG ] [ 0 ] );    
      
      // Create the NB term.
      
      double nbTerm = para.NB_BASE;
      
      nbTerm += initialThreshold * para.THRESHOLD_TO_NB;
      nbTerm += initialVoltage * para.VOLTAGE_TO_NB;
      
      sta_p [ State.NB ] = new State [ ] {
              new State ( 0, nbTerm, State.NB , 0 )
      };
      
      state.add ( sta_p [ State.NB ] [ 0 ] );
      
      // Create the NJ_SPIKE terms.

      sta_p [ State.NJ_SPIKE ] =
          new State [ para.allDecays [ State.NJ_SPIKE ].length ];
      
      for ( int  i = 0; i < sta_p [ State.NJ_SPIKE ].length; i++ )
      {
        double njTerm = 
          initialSpikeCurrents [ i ] * para.SPIKE_CURRENT_TO_SPIKE_NJ [ i ];
        
        sta_p [ State.NJ_SPIKE ] [ i ] =
                new State ( 0 , njTerm , State.NJ_SPIKE , i );
        
        state.add ( sta_p [ State.NJ_SPIKE ] [ i ] );
      }
      
      // Create the NJ_SYNAPSE terms.
      
      sta_p [ State.NJ_SYNAPSE ] = 
          new State [ para.allDecays [ State.NJ_SYNAPSE ].length ];
      
      for ( int  i = 0; i < sta_p [ State.NJ_SYNAPSE ].length; i++ )
      {
        sta_p [ State.NJ_SYNAPSE ] [ i ] = 
                new State ( 0 , 0 , State.NJ_SYNAPSE , i );
        
        state.add ( sta_p [ State.NJ_SYNAPSE ] [ i ] );
      }

      // Schedule the first fire event, if one should exist.
      
      // Use a modified Newton-Raphson (NR) root-finding iteration method to 
      // find the first zero of the thresholdDiff function.

      double  deltaT, thresholdDiff;
      
      double  deriv = safeDerivative ( 0.0 );
      
      // Ever-increasing estimate of the absolute time of the next fire
      // event.
      
      double  nextTime = ( -( initialVoltage - initialThreshold ) / deriv );
      
      deltaT = nextTime;
            
      while ( deriv > 0 && nextTime < 1.0
        && !( deltaT > 0 ? deltaT < tEPS : -deltaT < tEPS ) ) 
      {
        // Safe derivative of the thresholdDiff function at time = nextTime.
          
        deriv = safeDerivative ( nextTime );
        
        // Value of the thresholdDiff function at time = nextTime.
        
        thresholdDiff = thresholdDiff ( nextTime );
        
        // As per the NR algorithm.
        
        deltaT = ( -thresholdDiff / deriv );
        
        nextTime += deltaT; 
      }

      // If the derivative is decreasing or if the algorithm has gone far
      // enough, conclude that this neuron will not fire.
      
      if ( deriv < 0 || nextTime > 1.0 )
      {
        fire = false;
        
        nextTime = -1.0;
      }
      
      // If deltaT is within the proper precision, schedule an actual fire 
      // event.

      else if ( ( deltaT > 0 ? deltaT < tEPS : -deltaT < tEPS ) )
      {
        fire = true;
      }
      else
      {
        // This shouldn't happen.
      }

      if ( nextTime < 0 )
      {
      // A negative value for nextTime means that this neuron will not fire.
      }
      else
      {
        final Slot<FireEvent>  fireEventSlot = net.getFireEventSlot ( );
        
        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 ( )
    ////////////////////////////////////////////////////////////////////////
    {      
    /*
     * I.  If this is a real fire event:
     *      A.  Update spike time dependent plasticity variables. 
     *      B.  Update the state variables to time at the end of the
     *          refractory period.
     *      C.  Reset the membrane voltage, membrane voltage threshold, and 
     *          spike-induced currents according to the reset rules of the
     *          Mihalas-Niebur model.
     * II. Schedule the next fire event.
     */
      // The time to start looking for the next fire event.
        
      double  baselineTime;
      
      // The voltage and threshold at baselineTime.
      
      double  nowVoltage, nowThreshold;
      
      if ( fire ) // Only proceed if this is an actual fire event.
      {
        
        // Update STDP variables.
          
        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;
        
        // Back up the current threshold.
        
        final double  thresholdBackup
          = membraneVoltage ( timeOfNextFire ) 
          - thresholdDiff ( timeOfNextFire );
        
        // Update the G term to the time at the end of the refractory period.
                        
        sta_p [ State.NG ] [ 0 ].value *= Math.exp (
          -( timeOfNextFire - sta_p [ State.NG ] [ 0 ].time + para.ABSREF )
          * para.allDecays [ State.NG ] [ 0 ]);
        
        sta_p [ State.NG ] [ 0 ].time = timeOfNextFire + para.ABSREF;
        
        // Update the B term to the time at the end of the refractory period.
        
        sta_p [ State.NB ] [ 0 ].value *= Math.exp (
          -( timeOfNextFire - sta_p [ State.NB ] [ 0 ].time + para.ABSREF )
          * para.allDecays [ State.NB ] [ 0 ]);
        
        sta_p [ State.NB ] [ 0 ].time = timeOfNextFire + para.ABSREF;

        // Reset the membrane voltage and voltage threshold.
        
        // Calculate what the voltage and threshold will be
        // at the end of the refractory period.
        
        final double
          futureVoltage = membraneVoltage ( timeOfNextFire + para.ABSREF );
                
        final double 
          futureThresholdGap = thresholdDiff (timeOfNextFire + para.ABSREF);
        
        final double futureThreshold = futureVoltage - futureThresholdGap;
        
        // Calculate the reset voltage and threshold, according 
        // to the reset rules of the model.
        
        final double resetVoltage = para.VRESET;
        
        final double resetThreshold = 
                Math.max(para.RRESET, thresholdBackup + para.THRESHOLDADD);
        
        // Calculate the necessary offsets.
        
        final double  voltageOffset = resetVoltage - futureVoltage;
                
        final double  thresholdOffset = resetThreshold - futureThreshold;
        
        // Incorporate the offsets into the state variables.

        sta_p [ State.NG ] [ 0 ].value +=
            voltageOffset * para.VOLTAGE_TO_NG;
        
        sta_p [ State.NB ] [ 0 ].value +=
            voltageOffset * para.VOLTAGE_TO_NB;
        
        sta_p [ State.NB ] [ 0 ].value +=
            thresholdOffset * para.THRESHOLD_TO_NB; 
        
        // Clamp the voltage and threshold to the reset values.

        clampVoltage = resetVoltage;
        
        clampThreshold = resetThreshold;
        
        // Reset the spike-induced currents.
        
        for (int  a = 0; a < para.allDecays [ State.NJ_SPIKE ].length; a++)
        {
          // Reset a current only if its reset rule is not the identity.
          
          if ( para.SPIKE_RATIO [ a ] != 1.0 || para.SPIKE_ADD [ a ] != 0 )
          {            
            // Update the current to the end of the refractory period.
            
            // If the current is inactive, initialize its state variable.
                  
            if ( sta_p [ State.NJ_SPIKE ] [ a ] == null )
            {              
              sta_p [ State.NJ_SPIKE ] [ a ] = new State (
                timeOfNextFire + para.ABSREF,
                0,
                State.NJ_SPIKE,
                a);
              
              state.add ( sta_p [ State.NJ_SPIKE ] [ a ]);
            }
            else
            {
              sta_p [ State.NJ_SPIKE ] [ a ].value *= Math.exp (
                -( timeOfNextFire - 
                    sta_p [ State.NJ_SPIKE ] [ a ].time + para.ABSREF )
                * para.allDecays [ State.NJ_SPIKE ] [ a ] );
              
              sta_p [ State.NJ_SPIKE ] [ a ].time =
                  timeOfNextFire + para.ABSREF;
            }

            // Calculate what this current will be
            // at the end of the refractory period.
            
            final double  futureCurrent = 
                sta_p [ State.NJ_SPIKE ] [ a ].value
                    / para.SPIKE_CURRENT_TO_SPIKE_NJ [ a ];
            
            // Calculate the reset current, according to the reset rule 
            // of the model.
                        
            final double  resetCurrent =  
                ( para.SPIKE_RATIO [ a ]* futureCurrent
                    * Math.exp ( para.ABSREF *
                  para.allDecays [ State.NJ_SPIKE ] [ a ] )
              + para.SPIKE_ADD [ a ] )
              * Math.exp ( -para.ABSREF *
                  para.allDecays [ State.NJ_SPIKE ] [ a ] );
            
            // Calculate the necessary offset.
            
            final double  currentOffset = resetCurrent - futureCurrent;
            
            // Incorporate the offset into the state variables.
            
            sta_p [ State.NG ] [ 0 ].value += 
                currentOffset * para.SPIKE_CURRENT_TO_NG [ a ];
            
            sta_p [ State.NB ] [ 0 ].value +=
                currentOffset * para.SPIKE_CURRENT_TO_NB [ a ];
            
            sta_p [ State.NJ_SPIKE ] [ a ].value += 
                currentOffset * para.SPIKE_CURRENT_TO_SPIKE_NJ [ a ];
          }
        }
        
        timeOfLastUpdate = timeOfNextFire + para.ABSREF;
        
        baselineTime  = timeOfNextFire + para.ABSREF;
        
        nowVoltage = resetVoltage;
        
        nowThreshold = resetThreshold;
      }
      else // If this fire event was only a preliminary prediction.
      {
        baselineTime = timeOfNextFire;
        
        nowVoltage = thresholdDiff(baselineTime);
        
        nowThreshold = 0;
      }
      
      // Schedule the next fire event, if one should exist.
      
      // Use a modified Newton-Raphson (NR) root-finding iteration method to 
      // find the first zero of the thresholdDiff function.
      
      double  deltaT, thresholdDiff;
      
      double  deriv = safeDerivative ( baselineTime );
      
      // Ever-increasing estimate of the time until the next fire event.
      
      double  nextTime = ( -( nowVoltage - nowThreshold ) ) / deriv;
      
      deltaT = nextTime;
      
      // Repeat the NR iteration until either the derivative turns negative
      // or the predicted fire time grows so far away that an input spike
      // is likely to arrive in the intervening time, rendering any further
      // calculations useless (see D'Haene, 2009).

      while ( deriv > 0 && nextTime < cost * tAvg ) 
      {
        // Safe derivative of the thresholdDiff function.
          
        deriv = safeDerivative ( baselineTime + nextTime );
        
        // Value of the thresholdDiff function.
        
        thresholdDiff = thresholdDiff ( baselineTime + nextTime );
        
        // As per the NR algorithm.
        
        deltaT = ( -thresholdDiff / deriv );
        
        // If deltaT is within the proper precision, schedule an actual 
        // fire event.
        
        if ( deltaT > 0 ? deltaT < tEPS : -deltaT < tEPS )
        {
          fire = true;
                   
          return nextTime + baselineTime - timeOfNextFire;
        }
        
        nextTime += deltaT; 
      }
      
      // If the derivative is decreasing or if the algorithm has gone far 
      // enough, conclude that this neuron will not fire.

      if ( deriv < 0 || nextTime > 1.0 )
      {
        fire = false;
        
        return -1.0;
      }
      
      // Otherwise, schedule a preliminary estimate fire event.
      
      else
      {
        fire = false;
        
        return  nextTime + baselineTime - timeOfNextFire;
      }
    }

    @Override
    public double  updateInput (
      final double   time,
      final Synapse  input )
    ////////////////////////////////////////////////////////////////////////
    {
    /*
     * I.   Update spike time dependent plasticity variables. 
     * II.  Update the appropriate state variables to the current time.
     * III. Add the current from the input synapse into the appropriate
     *      state variables.
     * IV.  Schedule the next fire event.
     */
        
      // Update STDP variables.
                
      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 ) ) );
          }
        }
        
        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 );
      }

      // Update the running average of mean time intervals in between spikes.
      
      tAvg = tAvg * 0.8 + ( time - lastInputTime ) *0.2;
      
      // Store the time of this spike.
      
      lastInputTime = time;

      // Update the G term to the current time.
      
      sta_p [ State.NG ] [ 0 ].value *= Math.exp (
        -( time - sta_p [ State.NG ] [ 0 ].time ) 
        * para.allDecays [ State.NG ] [ 0 ]);

      sta_p [ State.NG ] [ 0 ].time = time;
      
      // Update the B term to the current time.
      
      sta_p [ State.NB ] [ 0 ].value *= Math.exp (
        -( time - sta_p [ State.NB ] [ 0 ].time )
        * para.allDecays [ State.NB ] [ 0 ]);
      
      sta_p [ State.NB ] [ 0 ].time = time;

      // Update the appropriate current to the current time.
      
      // If the appropriate current's channel is inactive, initialize its 
      // state variable.
      
      int channel = input.getType ( );
      
      if ( sta_p [ State.NJ_SYNAPSE ] [ channel ] == null ) 
      {   
        sta_p [ State.NJ_SYNAPSE ] [ channel ] =
            new State ( time, 0, State.NJ_SYNAPSE, channel);
        
        state.add ( sta_p [ State.NJ_SYNAPSE ] [ channel ] );
      }
      else
      {
        sta_p [ State.NJ_SYNAPSE ] [ channel ].value *= 
                Math.exp ( -( time - 
                    sta_p [ State.NJ_SYNAPSE ] [ channel ].time )
                        * para.allDecays [ State.NJ_SYNAPSE ] [ channel ] );
        
        sta_p [ State.NJ_SYNAPSE ] [ channel ].time = time;
      }

      // Add the input spike to the G term state variable.
      
      sta_p [ State.NG ] [ 0 ].value += 
          input.getWeight ( ) * ( histTable.get ( input ).weight )
        * para.SYNAPSE_CURRENT_TO_NG [ channel ];
      
      // Add the input spike to the B term state variable.
      
      sta_p [ State.NB ] [ 0 ].value += 
          input.getWeight ( ) * ( histTable.get ( input ).weight )
        * para.SYNAPSE_CURRENT_TO_NB [ channel ];
      
      // Add the input spike to the appropriate current's state variable. 
      
      sta_p [ State.NJ_SYNAPSE ] [ 0 ].value
      += input.getWeight ( ) 
          * histTable.get ( input ).weight
          * para.SYNAPSE_CURRENT_TO_SYNAPSE_NJ [ channel ];
      
      // The time to start looking for the next fire event.

      double baselineTime;

      // The voltage and threshold at baselineTime.
            
      double nowVoltage = membraneVoltage(time);
      
      double nowThreshold = nowVoltage - thresholdDiff(time);
      
      // If this neuron is still within a refractory period, offset the 
      // membrane voltage and threshold such that by the end of the 
      // refractory period, the voltage will be equal to the clamp voltage.
      
      if ( time < timeOfLastUpdate )
      {

        // Calculate what the voltage and threshold will be
        // at the end of the refractory period.
        
        double futureVoltage = membraneVoltage ( timeOfLastUpdate );
        
        double futureThresholdGap = thresholdDiff  ( timeOfLastUpdate );
        
        double futureThreshold = futureVoltage - futureThresholdGap;
        
        // Calculate the necessary offsets.
        
        double voltageOffset = clampVoltage - futureVoltage;
        
        double thresholdOffset = clampThreshold - futureThreshold;
        
        // Incorporate the offsets into the state variables.
        
        double ngDecay = Math.exp ( ( timeOfLastUpdate - time ) 
            * para.allDecays [ State.NG ] [ 0 ] );
        
        double nbDecay = Math.exp ( ( timeOfLastUpdate - time )
            * para.allDecays [ State.NB ] [ 0 ] );
        
        sta_p [ State.NG ] [ 0 ].value += para.VOLTAGE_TO_NG * voltageOffset 
            * ngDecay;
        
        sta_p [ State.NB ] [ 0 ].value += para.VOLTAGE_TO_NB * voltageOffset
            * nbDecay;
        
        sta_p [ State.NB ] [ 0 ].value += para.THRESHOLD_TO_NB
            * thresholdOffset * nbDecay;        
        
        nowVoltage = clampVoltage;
        
        nowThreshold = clampThreshold;
        
        baselineTime = timeOfLastUpdate;
      }
      else
      {
        clampVoltage = Double.MAX_VALUE;
        
        clampThreshold = Double.MAX_VALUE;
        
        timeOfLastUpdate = time;
        
        baselineTime = time;
      }
      
      // Schedule the next fire event, if one should exist.
      
      // Use a modified Newton-Raphson (NR) root-finding iteration method to 
      // find the first zero of the thresholdDiff function.

      double  deltaT, thresholdDiff;
      
      double  deriv = safeDerivative ( baselineTime );
      
      // Ever-increasing estimate of the time until the next fire event.

      double nextTime = ( -( nowVoltage - nowThreshold ) / deriv );
      
      deltaT = nextTime;
      
      // Repeat the NR iteration until either the derivative turns negative
      // or the predicted fire time grows so far away that an input spike
      // is likely to arrive in the intervening time, rendering any further
      // calculations useless (see D'Haene, 2009).
              
      while ( deriv > 0 && nextTime < cost * tAvg ) 
      {

          // Safe derivative of the thresholdDiff function.

          deriv = safeDerivative ( baselineTime + nextTime );

          // Value of the thresholdDiff function.

          thresholdDiff = thresholdDiff ( baselineTime + nextTime );

          // As per the NR algorithm.

          deltaT = ( -thresholdDiff / deriv );

          // If deltaT is within the proper precision, schedule an actual
          // fire event.

          if ( deltaT > 0 ? deltaT < tEPS : -deltaT < tEPS )
          {
              fire = true;

              return ( time > timeOfLastUpdate
                      ? nextTime : nextTime + timeOfLastUpdate - time );
          }

          nextTime += deltaT; 
      }

      // If the derivative is decreasing or if the algorithm has gone far
      // enough, conclude that this neuron will not fire.

      if( deriv < 0 || nextTime > 1.0 )
      {        
          fire = false;

          return -1.0;
      }

      // Otherwise, schedule a preliminary estimate fire event.

      else
      { 
          fire = false;

          return ( time > timeOfLastUpdate
                  ? nextTime : nextTime + timeOfLastUpdate - time );
      }
    }


    ////////////////////////////////////////////////////////////////////////
    // accessor methods
    ////////////////////////////////////////////////////////////////////////

    public double  getSensoryWeight ( )
    ////////////////////////////////////////////////////////////////////////
    {
      throw new RuntimeException ( "This neuron type doesn't use the "
        + "Sensory Weight Functions!" );
    }

    public double  getTimeOfLastUpdate ( )
    ////////////////////////////////////////////////////////////////////////
    {
      return this.timeOfLastUpdate;
    }


    ////////////////////////////////////////////////////////////////////////
    // 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+"Current:"+"\n";
      
      return tmp;
    }

    ////////////////////////////////////////////////////////////////////////
    // miscellaneous methods
    ////////////////////////////////////////////////////////////////////////

    
    /***********************************************************************
     * Compute the membrane voltage.
     * 
     * This method transforms the thresholdDiff (Mihalas and Niebur, 2009,
     *  Equation 3.5) into the voltage ( V(t) in Equation 3.2)
     * 
     * @param t
     *            absolute time t
     * @return membrane voltage at time t.s
     **********************************************************************/
    public double  membraneVoltage ( final double  t )
    ////////////////////////////////////////////////////////////////////////
    {

      double  V = 0;
      
      double  constant = 0;
      
      Iterator<State>  iter = state.iterator ( );
      
      State  first = iter.next ( );
      
      State  second;
      
      double  firstV, secondV;

      if ( first.type == State.NB ) // Ignore the NB term.
      {
        if ( iter.hasNext ( ) )
        {
          first = iter.next ( );
        }
        else
        {
          return V + constant + para.MEMBRANE_VOLTAGE_BASE + para.VREST;
        }
      }

      if ( first.time == t ) // Ignore no time change term. 
        // TODO: what does that mean?
      {
        if ( first.type == State.NB )
        {
          constant += first.value / para.VOLTAGE_TO_NG;
        }
        else if (first.type == State.NJ_SPIKE)
        {
          constant += first.value 
              / para.VOLTAGE_TO_SPIKE_NJ [ first.index ];
        }
        else if (first.type == State.NJ_SYNAPSE)
        {
          constant += first.value 
              / para.VOLTAGE_TO_SYNAPSE_NJ [ first.index ];
        }
        
        if ( iter.hasNext ( ) )
        {
          first = iter.next ( );
        }
        else
        {
          return V + constant + para.MEMBRANE_VOLTAGE_BASE + para.VREST;
        }
      }

      if ( first.type == State.NB ) // Ignore the NB term.
      {
        if ( iter.hasNext ( ) )
        {
          first = iter.next ( );
        }
        else
        {
          return V + constant + para.MEMBRANE_VOLTAGE_BASE + para.VREST;
        }
      }

      if ( first.type == State.NG )
      {
        firstV = first.value / para.VOLTAGE_TO_NG;
      }
      else if (first.type == State.NJ_SPIKE)
      {        
        firstV = first.value / para.VOLTAGE_TO_SPIKE_NJ [ first.index ];
      }
      else
      {        
        firstV = first.value / para.VOLTAGE_TO_SYNAPSE_NJ [ first.index ];
      }

      V = firstV;
      
      while ( iter.hasNext ( ) )
      {
        second = iter.next ( );
        
        secondV = 0.0;

        if ( second.time == t && second.type != State.NB )
        {
          if ( second.type == State.NG )
          {
            constant += second.value / para.VOLTAGE_TO_NG;
          }
          else if (second.type == State.NJ_SPIKE)
          {
            constant += second.value 
                / para.VOLTAGE_TO_SPIKE_NJ [ second.index ];
          }
          else
          {
            constant += second.value
                / para.VOLTAGE_TO_SYNAPSE_NJ [ second.index ];
          }
          
          if ( iter.hasNext ( ) )
          {
            second = iter.next ( );
          }
          else
          {
            break;
          }
        }

        if ( second.type == State.NB )
        {
          if ( iter.hasNext ( ) )
          {
            second = iter.next ( );
          }
          else
          {
            break;
          }
        }

        if ( second.type == State.NG )
        {
          secondV = second.value / para.VOLTAGE_TO_NG;
        }
        else if (second.type == State.NJ_SPIKE)
        {
          secondV = second.value
              / para.VOLTAGE_TO_SPIKE_NJ [ second.index ];
        }
        else 
        {
          secondV = second.value 
              / para.VOLTAGE_TO_SYNAPSE_NJ [ second.index ];
        }

        V = Math.exp ( -( ( para.allDecays [ first.type ] [ first.index ]
          - para.allDecays [ second.type ] [ second.index ] ) * t
          - ( para.allDecays [ first.type ] [ first.index ] * first.time
          - para.allDecays [ second.type ] [ second.index ] * second.time )
          ) ) * V + secondV;
     
        first = second;
      }

      V = Math.exp ( -( para.allDecays [ first.type ] [ first.index ] * t
        - para.allDecays [ first.type ] [ first.index ] * first.time )) * V;
            
      return V + para.MEMBRANE_VOLTAGE_BASE + para.VREST + constant;
    }
    
    /***********************************************************************
     * Compute the difference between the membrane voltage and the membrane
     * voltage threshold. The <code>MNNeuron</code> class is optimized to
     * perform this function quickly, at the expense of
     * <code>membraneVoltage()</code>.
     * 
     * This method is based on Equation 3.5 in Mihalas and Niebur, 2009.
     * 
     * @param t
     *            absolute time t
     * @return voltage - threshold at time t
     **********************************************************************/
    public double  thresholdDiff ( double  t )
    ////////////////////////////////////////////////////////////////////////
    {
      double  V = 0, constant = 0.0;
      
      Iterator<State>  iter = state.iterator ( );
      
      State  first = iter.next ( );
      State  second;
      
      double  firstV, secondV;
      
      // If the state is a current channel and its value is too small, delete
      // it.

      if ( Math.abs ( first.value ) < rEPS
                  && ( int ) first.type != State.NG 
                  && ( int ) first.type != State.NB ) 
      {        
        iter.remove ( );
        
        sta_p [ first.type ] [ first.index ] = null;
      }
      
      // ignore no time change term
      
      if ( first.time == t )
      {
        constant += first.value;
        
        if ( iter.hasNext ( ) )
        {
          first = iter.next ( );
        }
        else
        {
          return V + constant + para.THRESHOLD_DIFF_BASE;
        }
      }
      
      firstV = first.value;
      
      V = firstV;
      
      while ( iter.hasNext ( ) )
      {
        second = iter.next ( );
        
        // If the state is a current channel and its value is too small,
        // delete it.
        
        if ( Math.abs ( second.value ) < rEPS 
                  && second.type != State.NG 
                  && second.type != State.NB ) 
        {
          iter.remove ( );
          
          sta_p [ second.type ] [ second.index ] = 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 [ first.type ] [ first.index ]
          - para.allDecays [ second.type ] [ second.index ] ) * t
          - ( para.allDecays [ first.type ] [ first.index ] * first.time
          - para.allDecays [ second.type ] [ second.index ] * second.time )
          ) ) * V + secondV;
        
        first = second;
      }

      V = Math.exp ( -( para.allDecays [ first.type ] [ first.index ] * t
        - para.allDecays [ first.type ] [ first.index ] * first.time )) * V;
            
      return V + para.THRESHOLD_DIFF_BASE + constant;
    }

    /***********************************************************************
    * Computes the "safe derivative" (see D'Haene, 2009) of the 
    * thresholdDiff function. 
    * 
    * @param t
    *   absolute time t
    ***********************************************************************/
    public double  safeDerivative ( 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();

        if ( tmpState.value < 0 )
        {
          tauSafe = para.allDecays [ tmpState.type ] [ tmpState.index];
          
          break;
        }
      }

      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 [ first.type ] [ first.index ],
            tauSafe );
          
          constant += -first.value * maxDecay;
        }
        else
        {
          constant += 
              -first.value * para.allDecays [ first.type ] [ first.index ];
        }

        if ( iter.hasNext ( ) )
        {
          first = iter.next ( );
        }
        else
        {
          return V + constant;
        }
      }

      if ( first.value > 0 )
      {
        final double  maxDecay = Math.min (
          para.allDecays [ first.type ] [ first.index ],
          tauSafe );
        
        firstV = -first.value * maxDecay;
      }
      else
      {
        firstV = 
            -first.value * para.allDecays [ first.type ] [ first.index ];
      }

      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 [ second.type ] [ second.index ],
              tauSafe );
            
            constant += -second.value * maxDecay;
          }
          else
          {
            constant
              += -second.value
              * para.allDecays [ second.type ] [ second.index ];
          }

          if ( iter.hasNext ( ) )
          {
            second = iter.next ( );
          }
          else
          {
            break;
          }
        }

        if ( second.value > 0 )
        {
          final double  maxDecay = Math.min (
            para.allDecays [ second.type ] [ second.index ],
            tauSafe );
          
          secondV = -second.value * maxDecay;
        }
        else
        {
          secondV =
              -second.value * para.allDecays [ second.type] [ second.index];
        }

        V = Math.exp ( -( ( para.allDecays [ first.type ] [ first.index]
          - para.allDecays [ second.type ] [ second.index] ) * t
          - ( para.allDecays [ first.type ] [ first.index ] * first.time
          - para.allDecays [ second.type ] [ second.index] * second.time ) 
          ) ) * V + secondV;
        
        first = second;
      }

      V = Math.exp ( -( para.allDecays [first.type ] [ first.index ] * t
        - para.allDecays [ first.type ] [ first.index ] * first.time )) * V;
      
      return V + constant;
    }

    ////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////
    }