SpikingNeuron.cs

Go to the documentation of this file.

namespace SpikingNeuronNetwork.Lib.NeuronModels
{
    using Interfaces;
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;
    using Training;

    public abstract class SpikingNeuron : ISpikingNeuron
    {
        public static double Timestep = 0.0001;
        
        public static double Epsilon = 1e-10;
        
        public const int MaxNumOutputSpikes = 5;

        protected Dictionary<Synapse, double> _weights;
        
        protected SpikePriorityQueue SpikeQueue;
        
        private readonly Random _random = new Random();

        public int NeuronIndex { get; set; }

        public SpikingNeuronNetwork NeuronNetwork { get; set; }

        public bool Verbose { get; set; }

        public NeuronState State { get; protected set; }

        public SimulationMethod Method { get; set; }

        public abstract double PostSpikeState { get; }

        public abstract double SpikeThreshold { get; }

        public abstract double InitialState { get; }

        public abstract double RemainingTimeToSpike { get; }

        public abstract double MaxNumericalIntegrationIterations { get; }
        
        public abstract OperatingRegion OperatingRegion { get; }

        public double NextSpikeTime
        {
            get { return RemainingTimeToSpike + State.Time; }
        }

        public bool IsStandAloneNeuron
        {
            get { return NeuronNetwork == null; }
        }

        public int NumInputs
        {
            get
            {
                return IsStandAloneNeuron
                           ? Weights.Count(x => Math.Abs(x.Value) > Epsilon)
                           : Weights.Where(x => x.Key.InputNeuronIndex < NeuronNetwork.NumInputs)
                                    .Count(x => Math.Abs(x.Value) > Epsilon);
            }
        }

        public Dictionary<Synapse, double> Weights
        {
            get { return IsStandAloneNeuron ? _weights : NeuronNetwork.GetNeuronWeights(NeuronIndex); }
        }

        protected SpikingNeuron()
            : this(0)
        {
        }

        protected SpikingNeuron(SpikingNeuronNetwork network, int neuronIndex)
        {
            NeuronNetwork = network;
            NeuronIndex = neuronIndex;
            ResetState();
            Verbose = true;
            Method = SimulationMethod.EventDriven;
        }

        protected SpikingNeuron(IReadOnlyList<double> weights)
        {
            var numInputs = weights.Count;
            NeuronIndex = numInputs;
            _weights = new Dictionary<Synapse, double>(weights.Count);
            for (var k = 0; k < numInputs; k++)
            {
                _weights.Add(new Synapse(k, NeuronIndex), weights[k]);
            }
            ResetState();
            Verbose = true;
            Method = SimulationMethod.Numerical;
        }

        protected SpikingNeuron(int numInputs = 3, double weightIni = 1e6)
        {
            _weights = new Dictionary<Synapse, double>(numInputs);
            NeuronIndex = numInputs;
            for (var k = 0; k < numInputs; k++)
            {
                var newWeight = weightIni > 100 ? 2*_random.NextDouble() - 1 : weightIni;
                _weights.Add(new Synapse(k, NeuronIndex), newWeight);
            }
            ResetState();
            Verbose = true;
            Method = SimulationMethod.Numerical;
        }

        public abstract ISpikingNeuron Clone();

        public abstract List<Spike> RunThetaNeuron(List<Spike> inputSpikes, int maxNumOutputSpikes = MaxNumOutputSpikes);

        public abstract IEnumerable<Spike> UpdateStateVariableOnSpike(double weight);

        public abstract IEnumerable<Spike> AdvanceNeuronState(double newTime);

        public abstract List<double> CalculatePostSpikeDerivatives(NeuronFiringHistory neuronFiringHistory);

        public abstract Dictionary<Synapse, NeuronDerivativeParameters> CalculateOutputSpikeTimeDerivatives(NeuronFiringHistory neuronFiringHistory);
        
        public abstract double StateDerivative(double stateVariable);

        public List<Spike> RunThetaNeuronNumerically(int maxNumOutputSpikes = MaxNumOutputSpikes)
        {
            var outputSpikes = new List<Spike>();
            if (SpikeQueue.Count == 0)
            {
                return outputSpikes;
            }
            var statePast = State.StateVariable;
            var allInputSpikeProcessed = false;
            var currentSpike = SpikeQueue.Dequeue();

            //The maximum number of numerical integration iterations
            var intLimit = MaxNumericalIntegrationIterations;

            var doIntegration = true;
            while (doIntegration)
            {
                // Forward Euler Numeric Integration
                State.StateVariable = statePast + Timestep * StateDerivative(statePast);

                if (!allInputSpikeProcessed && (currentSpike.Time >= State.Time) && (currentSpike.Time < (State.Time + Timestep)))
                {
                    var synapse = new Synapse(currentSpike.NeuronIndex, NeuronIndex);
                    var spikeStats = new SpikeStats();
                    try 
                    {
                        spikeStats = new SpikeStats
                            {
                                PreSpikeState = new NeuronState {StateVariable = State.StateVariable, Time = currentSpike.Time},
                                Weight = Weights[synapse],
                                Synapse = synapse
                            };
                    }
                    catch (KeyNotFoundException)
                    {
                        Console.WriteLine("Synapse: " + synapse);
                        Console.WriteLine("Weights: " + GetWeightsString(Weights));
                        throw;
                    }
                    UpdateStateVariableOnSpike(spikeStats.Weight);
                    spikeStats.PostSpikeState = new NeuronState { StateVariable = State.StateVariable, Time = currentSpike.Time };
                    if (Verbose)
                    {
                        Console.WriteLine("*** Spike From Neuron " + currentSpike.NeuronIndex + " To Neuron " + NeuronIndex + " ***");
                        Console.WriteLine(spikeStats);
                    }
                    if (SpikeQueue.Count > 0)
                    {
                        currentSpike = SpikeQueue.Dequeue();
                    }
                    else
                    {
                        allInputSpikeProcessed = true;
                    }
                }

                // If all input spikes have been processed and phase is quiescent
                // then no further output spikes are possible
                if (allInputSpikeProcessed && (OperatingRegion != OperatingRegion.SpikeGeneration && OperatingRegion != OperatingRegion.TonicSpiking))
                {
                    break;
                }

                // Check for output spike
                if (State.StateVariable > SpikeThreshold && statePast <= SpikeThreshold)
                {
                    outputSpikes.Add(new Spike { NeuronIndex = NeuronIndex, Time = State.Time });
                    if (Verbose)
                    {
                        Console.WriteLine("---> " + outputSpikes.Last() + "\n");
                    }
                    State.StateVariable = PostSpikeState;
                    if (outputSpikes.Count == maxNumOutputSpikes)
                    {
                        return outputSpikes;
                    }
                }

                // Special Case of Refractory Firing
                if (State.StateVariable > PostSpikeState && statePast <= PostSpikeState)
                {
                    State.StateVariable = SpikeThreshold - (State.StateVariable + PostSpikeState);
                }

                // Prepare for Next Iteration
                statePast = State.StateVariable;
                State.Time += Timestep;
                if (State.Time <= intLimit) continue;
                if (!(OperatingRegion == OperatingRegion.SpikeGeneration ||
                        (OperatingRegion == OperatingRegion.TonicSpiking && (outputSpikes.Count == 0 || currentSpike.Time > outputSpikes.Last().Time))))
                {
                    doIntegration = false;
                }
            }

            return outputSpikes;
        }

        public List<double> CalculatePostSpikeDerivativesNumerical(NeuronFiringHistory neuronFiringHistory)
        {
            // Calculate all DthpDthpm1 terms (Equation 3.1.4 in NECO Paper, Although MIMO Version is Used Here)
            var postSpikeDerivativeProducts = new List<double>();
            SpikeStats previousSpikeStat = null;
            var initialState = State;
            const double offset = 0.001;
            foreach (var spikeStat in neuronFiringHistory.SpikeStats)
            {
                if (previousSpikeStat == null)
                {
                    previousSpikeStat = spikeStat;
                    continue;
                }

                var currentSpike = new Spike
                {
                    NeuronIndex = spikeStat.Synapse.InputNeuronIndex,
                    Time = spikeStat.PreSpikeState.Time
                };
                
                State = new NeuronState { StateVariable = previousSpikeStat.PostSpikeState.StateVariable, Time = previousSpikeStat.PostSpikeState.Time };
                State.StateVariable -= offset;
                SpikeStats spikeStatsGradDelta1;
                ProcessSpike(currentSpike, out spikeStatsGradDelta1);
                State = new NeuronState { StateVariable = previousSpikeStat.PostSpikeState.StateVariable, Time = previousSpikeStat.PostSpikeState.Time };
                State.StateVariable += offset;
                SpikeStats spikeStatsGradDelta2;
                ProcessSpike(currentSpike, out spikeStatsGradDelta2);

                postSpikeDerivativeProducts.Add((spikeStatsGradDelta2.PostSpikeState.StateVariable - spikeStatsGradDelta1.PostSpikeState.StateVariable) / (2 * offset));
                previousSpikeStat = spikeStat;
            }

            State = initialState;
            return postSpikeDerivativeProducts;
        }

        public Dictionary<Synapse, NeuronDerivativeParameters> CalculateOutputSpikeTimeDerivativesNumerical(NeuronFiringHistory neuronFiringHistory)
        {
            // Convert neuronFiringHistory to input spikes
            var inputSpikes = neuronFiringHistory.SpikeStats.Select(spikeStat => new Spike
            {
                NeuronIndex = spikeStat.Synapse.InputNeuronIndex,
                Time = spikeStat.PreSpikeState.Time
            }).ToList();

            if (neuronFiringHistory.OutputSpikes.Count > 1)
            {
                throw new ArgumentException("MIMO Training Not Implemented");
            }

            if (neuronFiringHistory.OutputSpikes.Count == 0)
            {
                throw new ArgumentException("Actual Output Spike Times are Missing");
            }

            var inputSynapses = Weights.Select(x => x.Key).ToList();
            var outputSpikeTimeDerivatives = inputSynapses.ToDictionary(synapse => synapse, synapse => new NeuronDerivativeParameters
                {
                    OutputSpikeTimeToInputSpikeTimeDerivative = 0, OutputSpikeTimeToWeightDerivative = 0
                });

            const double offset = 0.00001;
            var initialState = State;
            foreach (var spikeStat in neuronFiringHistory.SpikeStats)
            {
                // Modify weight
                var weight = Weights.First(x => x.Key.InputNeuronIndex == spikeStat.Synapse.InputNeuronIndex);
                
                SetSynapticWeight(spikeStat.Synapse, weight.Value - offset);
                ResetState();
                var outputSpikeTime1 = RunThetaNeuron(inputSpikes).First().Time;
                
                SetSynapticWeight(spikeStat.Synapse, weight.Value + offset);
                ResetState();
                var outputSpikeTime2 = RunThetaNeuron(inputSpikes).First().Time;

                SetSynapticWeight(spikeStat.Synapse, weight.Value);
                var outputToWeightDerivative = (outputSpikeTime2 - outputSpikeTime1) / (2 * offset);

                // Modify input spike time
                var currentSpike = inputSpikes.First(x => x.NeuronIndex == spikeStat.Synapse.InputNeuronIndex);

                currentSpike.Time -= offset;
                ResetState();
                var outputSpikeTime3 = RunThetaNeuron(inputSpikes).First().Time;

                currentSpike.Time += 2*offset;
                ResetState();
                var outputSpikeTime4 = RunThetaNeuron(inputSpikes).First().Time;

                var outputToInputSpikeTimeDerivative = (outputSpikeTime4 - outputSpikeTime3) / (2 * offset);

                var derivativeParameters = new NeuronDerivativeParameters
                    {
                        OutputSpikeTimeToWeightDerivative = outputToWeightDerivative,
                        OutputSpikeTimeToInputSpikeTimeDerivative = outputToInputSpikeTimeDerivative
                    };
                outputSpikeTimeDerivatives[spikeStat.Synapse] = derivativeParameters;
            }

            State = initialState;
            return outputSpikeTimeDerivatives;
        }

        public IEnumerable<Spike> ProcessSpike(Spike spike, out SpikeStats spikeStats)
        {
            // Protect against non-causality
            if (State.Time > spike.Time)
            {
                spikeStats = null;
                return new List<Spike>();
            }

            var outputSpikeTimes = AdvanceNeuronState(spike.Time).ToList();
            spikeStats = new SpikeStats
            {
                PreSpikeState = new NeuronState
                {
                    StateVariable = State.StateVariable,
                    Time = State.Time
                },
                Weight = Weights.First(x => x.Key.InputNeuronIndex == spike.NeuronIndex).Value,
                Synapse = new Synapse(spike.NeuronIndex, NeuronIndex)
            };

            var updateSpikes = UpdateStateVariableOnSpike(spikeStats.Weight).ToList();
            var firedOnSpike = updateSpikes.Count > 0;
            outputSpikeTimes.AddRange(updateSpikes);

            spikeStats.PostSpikeState = new NeuronState
            {
                StateVariable = State.StateVariable,
                Time = State.Time
            };

            if (Verbose)
            {
                Console.WriteLine("*** Spike From Neuron " + spike.NeuronIndex + " To Neuron " + NeuronIndex + " ***");
                Console.WriteLine(spikeStats);
                if (outputSpikeTimes.Count > 0)
                {
                    foreach (var outputSpikeTime in outputSpikeTimes)
                    {
                        Console.WriteLine("---> " + outputSpikeTime);
                    }
                }
                else
                {
                    Console.WriteLine("---> No Output Spikes Produced");
                }

                if (firedOnSpike)
                {
                    Console.WriteLine("---> At least one output spike fired directly from input spike");
                }

                Console.WriteLine();
            }

            return outputSpikeTimes;
        }

        public void ResetState()
        {
            State = new NeuronState
                {
                    StateVariable = InitialState,
                    Time = 0
                };
        }

        public void SetSynapticWeight(Synapse synapse, double newValue)
        {
            if (IsStandAloneNeuron)
            {
                if (_weights.ContainsKey(synapse))
                {
                    _weights[synapse] = newValue;
                }
                else
                {
                    _weights.Add(synapse, newValue);    
                }                
            }
            else
            {
                NeuronNetwork.SetSynapticWeight(synapse, newValue);
            }
        }

        protected static string GetWeightsString(Dictionary<Synapse, double> weights)
        {
            var weightsStringBuilder = new StringBuilder();
            weights.Keys.ToList().ForEach(x => weightsStringBuilder.AppendLine("\t" + weights[x] + "(" + x.InputNeuronIndex + "->" + x.OutputNeuronIndex + ")"));
            return weightsStringBuilder.ToString();
        }        
    }
}