ThetaNeuron.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 class ThetaNeuron : SpikingNeuron
    {
        public const double DefaultIo = -0.001;

        public static double Alpha = 1;

        public static double Offset = 0.0001;

        public double Io { get; set; }

        #region Readonly Properties

        public override double PostSpikeState
        {
            get { return -Math.PI; }
        }

        public override double SpikeThreshold
        {
            get { return Math.PI; }
        }

        public override double InitialState
        {
            get { return (Io < 0) ? PositiveFixedPoint + Offset : 0; }
        }

        public override double RemainingTimeToSpike
        {
            get
            {
                if (Io > 0)
                {
                    return (1 / Beta) * ((SpikeThreshold / 2) - Math.Atan(C1));
                }
                if (State.StateVariable > PositiveFixedPoint && State.StateVariable < SpikeThreshold)
                {
                    return (1 / Beta) * Atanh(1 / C1);
                }
                return Double.PositiveInfinity;
            }
        }

        public override OperatingRegion OperatingRegion
        {
            get
            {
                if (Io > 0)
                {
                    return OperatingRegion.TonicSpiking;
                }
                if (State.StateVariable < PositiveFixedPoint && State.StateVariable > -PositiveFixedPoint)
                {
                    return OperatingRegion.Quiescent;
                }
                if (State.StateVariable >= PositiveFixedPoint)
                {
                    return OperatingRegion.SpikeGeneration;
                }
                return OperatingRegion.Refractory;
            }
        }

        public override double MaxNumericalIntegrationIterations
        {
            get
            {
                double simulationLimit;
                // Integrate long enough to include effects of all input spikes
                if (SpikeQueue.Count > 0)
                {
                    simulationLimit = (SpikeQueue.Last().Time + 1.5*(Io < 0 ? BaselineFiringTime : RemainingTimeToSpike));
                }
                else
                {
                    simulationLimit = 1.1*BaselineFiringTime;
                }
                return simulationLimit;
            }
        }

        public double Period
        {
            get { return Io > 0 ? SpikeThreshold/Beta : Double.NaN; }
        }

        public double Beta
        {
            get { return Math.Sqrt(Math.Abs(Alpha*Io)); }
        }

        public double PositiveFixedPoint 
        {
            get { return (Io < 0) ? Math.Acos((1 + Alpha*Io)/(1 - Alpha*Io)) : Double.NaN; }
        }

        public double BaselineFiringTime
        {
            get { return (1/Beta)*Atanh(Beta/Math.Tan((InitialState)/2)); }
        }

        private double C1
        {
            get { return Math.Tan(State.StateVariable/2)/Beta; }
        }

        #endregion

        private static double Atanh(double x)
        {
            return (0.5 * Math.Log((1.0 + x) / (1.0 - x)));
        }

        public ThetaNeuron()
            : this(0)
        {
        }

        public ThetaNeuron(SpikingNeuronNetwork network, int neuronIndex, double newIo = DefaultIo)
            : base(network, neuronIndex)
        {
            Io = newIo;
            ResetState();
        }

        public ThetaNeuron(IReadOnlyList<double> weights, double newIo = DefaultIo)
            : base(weights)
        {
            Io = newIo;
            ResetState();
        }

        public ThetaNeuron(int numInputs = 3, double weightIni = 1e6, double newIo = DefaultIo)
            : base(numInputs, weightIni)
        {
            Io = newIo;
            ResetState();
        }

        public override double StateDerivative(double phase)
        {
            return ((1 - Math.Cos(phase)) + Alpha*Io*(1 + Math.Cos(phase)));
        }

        public override ISpikingNeuron Clone()
        {
            return new ThetaNeuron(Weights.Select(x => x.Value).ToList())
                {
                    NeuronIndex = NeuronIndex,
                    NeuronNetwork = null,
                    State = State,
                    Method = Method,
                    Io = Io,
                    Verbose = Verbose,
                };
        }

        public override string ToString()
        {
            var neuronStringBuilder = new StringBuilder();
            neuronStringBuilder.AppendLine("\nNeuron Properties: ");
            neuronStringBuilder.AppendLine("\tBeta: " + String.Format("{0:0.0000}", Beta));
            neuronStringBuilder.AppendLine("\tIo: " + Io);
            neuronStringBuilder.AppendLine("\tPositive Fixed Point: " + String.Format("{0:0.0000}", PositiveFixedPoint));
            neuronStringBuilder.AppendLine("\tPhase: " + String.Format("{0:0.0000}", State.StateVariable));
            neuronStringBuilder.AppendLine("\tTime: " + String.Format("{0:0.0000}", State.Time));
            neuronStringBuilder.AppendLine("\tNumber of Inputs: " + NumInputs);
            neuronStringBuilder.AppendLine("\tWeights: ");
            neuronStringBuilder.Append("\t\t");
            neuronStringBuilder.AppendLine(GetWeightsString(Weights));
            neuronStringBuilder.AppendLine();
            neuronStringBuilder.AppendLine();
            return neuronStringBuilder.ToString();
        }   
   
        public override List<Spike> RunThetaNeuron(List<Spike> inputSpikes, int maxNumOutputSpikes = MaxNumOutputSpikes)
        {
            SpikeQueue = SpikePriorityQueue.CreateSpikeQueue(inputSpikes, State.Time);
            var outputSpikes = new List<Spike>();
            if (SpikeQueue.Count == 0)
            {
                return outputSpikes;
            }

            if (Verbose) 
            {
                Console.WriteLine("Input Spike Times: ");
                Console.Write(SpikeSet.GetSpikeTimesString(inputSpikes));
                Console.WriteLine("Baseline Firing Time: " + BaselineFiringTime);
                Console.WriteLine("Simulation Method: " + Method);
                Console.WriteLine();
            }

            switch (Method)
            {
                case SimulationMethod.EventDriven:
                    {
                        while (SpikeQueue.Count > 0)
                        {
                            SpikeStats spikeStats;
                            var newOutputSpikeTimes = ProcessSpike(SpikeQueue.Dequeue(), out spikeStats);
                            if (newOutputSpikeTimes != null)
                            {
                                outputSpikes.AddRange(newOutputSpikeTimes);
                            }

                            if (outputSpikes.Count>=maxNumOutputSpikes) 
                            {
                                return outputSpikes;
                            }
                        }

                        //Calculate thm later to save time, unless needed to spike checking

                        //If Io<0 and phase is between the positive fixed point and pi, calculate the final ts
                        //with or without previous input spikes
                        // If Io>0, produce an addition output spike if an input spike has occured since the 
                        //      last output spike or if no output spikes have been produced yet
                        if ((Io < 0 && State.StateVariable > PositiveFixedPoint) || (Io > 0 && (outputSpikes.Count == 0 || State.Time > outputSpikes.Last().Time)))
                        {
                            outputSpikes.Add(new Spike { NeuronIndex = NeuronIndex, Time = NextSpikeTime });
                            if (Verbose) 
                            {
                                Console.WriteLine("---> " + outputSpikes.Last());
                            }                                                        
                        }
                    }
                    break;

                case SimulationMethod.Numerical:
                    {
                        outputSpikes = RunThetaNeuronNumerically(maxNumOutputSpikes);
                    }
                    break;
            }

            if (Verbose)
            {
                Console.WriteLine();
            }

            return outputSpikes;
        }

        public override List<double> CalculatePostSpikeDerivatives(NeuronFiringHistory neuronFiringHistory)
        {
            var postSpikeDerivativeProducts = new List<double>();
            SpikeStats previousSpikeStat = null;
            foreach (var spikeStat in neuronFiringHistory.SpikeStats)
            {
                if (previousSpikeStat == null)
                {
                    previousSpikeStat = spikeStat;
                    continue;
                }

                postSpikeDerivativeProducts.Add(((1 + Math.Cos(spikeStat.PostSpikeState.StateVariable)) *
                    (Math.Pow(Math.Tan(spikeStat.PreSpikeState.StateVariable / 2), 2) + Alpha * Io)) /
                    StateDerivative(previousSpikeStat.PostSpikeState.StateVariable));
                previousSpikeStat = spikeStat;
            }

            return postSpikeDerivativeProducts;
        }

        public override Dictionary<Synapse, NeuronDerivativeParameters> CalculateOutputSpikeTimeDerivatives(NeuronFiringHistory neuronFiringHistory)
        {
            var nHatPostSpikeState = Double.NaN;
            var nHatPreSpikeState = Double.NaN;
            var nHatWeight = Double.NaN;
            var spikesProcessed = 0;
            var inputSynapses = Weights.Select(x => x.Key).ToList();

            var outputSpikeTimeDerivatives = inputSynapses.ToDictionary(synapse => synapse, synapse => new NeuronDerivativeParameters
                {
                    OutputSpikeTimeToInputSpikeTimeDerivative = 0, OutputSpikeTimeToWeightDerivative = 0
                });

            if (neuronFiringHistory.OutputSpikes.Count == 0)
            {
                return outputSpikeTimeDerivatives;
            }
            
            var actualFiringTime = neuronFiringHistory.OutputSpikes.First().Time;
            foreach (var spikeStat in neuronFiringHistory.SpikeStats.TakeWhile(spikeStat => spikeStat.PreSpikeState.Time < actualFiringTime))
            {
                nHatPreSpikeState = spikeStat.PreSpikeState.StateVariable;
                nHatPostSpikeState = spikeStat.PostSpikeState.StateVariable;
                nHatWeight = spikeStat.Weight;
                spikesProcessed++;
            }
            var nHat = spikesProcessed - 1;

            // No input spikes occur before the output spike, so little adjustments to weights will have no effect
            if (spikesProcessed == 0)
            {
                return outputSpikeTimeDerivatives;
            }

            // Simple Version when there is only one input spike
            if (neuronFiringHistory.SpikeStats.Count == 1)
            {
                var temp = (Math.Pow(Math.Tan(nHatPreSpikeState / 2), 2) + Alpha * Io);
                outputSpikeTimeDerivatives[neuronFiringHistory.SpikeStats.First().Synapse] = 
                    new NeuronDerivativeParameters
                    {
                        OutputSpikeTimeToInputSpikeTimeDerivative = temp / (temp + Math.Pow(Alpha * nHatWeight, 2) + 2 * Alpha * nHatWeight + Math.Tan(nHatPreSpikeState / 2)),
                        OutputSpikeTimeToWeightDerivative = -Alpha / (Math.Pow(Math.Tan(nHatPostSpikeState / 2), 2) + Alpha * Io)
                    };
                return outputSpikeTimeDerivatives;
            }

            // Multi-Synapse Version
            var nHatPostSpikeStateDerivative = StateDerivative(nHatPostSpikeState);
            var postSpikeDerivativeProducts = CalculatePostSpikeDerivatives(neuronFiringHistory);
            var n = 0;
            foreach (var spikeStat in neuronFiringHistory.SpikeStats)
            {
                var prod = 1.0;
                for (var j = n; j < nHat; j++)
                {
                    prod *= postSpikeDerivativeProducts[j];
                }
                var gamma = -spikeStat.Weight * (Alpha * spikeStat.Weight + 2 * Math.Tan(spikeStat.PreSpikeState.StateVariable / 2));
                var weightDerivative = ((-Alpha * (1 + Math.Cos(spikeStat.PostSpikeState.StateVariable))) / nHatPostSpikeStateDerivative) * prod;
                outputSpikeTimeDerivatives[spikeStat.Synapse]=
                    new NeuronDerivativeParameters
                    {
                        OutputSpikeTimeToWeightDerivative = weightDerivative,
                        OutputSpikeTimeToInputSpikeTimeDerivative = gamma * weightDerivative
                    };
                n++;
            }

            return outputSpikeTimeDerivatives;
        }        
        
        public override IEnumerable<Spike> UpdateStateVariableOnSpike(double weight)
        {
            // The following will approach pi, but never exceed it as Phase->pi and weight->Inf
            // So in effect, output spikes can never directly result from an input spike
            State.StateVariable = 2 * Math.Atan(Alpha * weight + Math.Tan(State.StateVariable / 2));
            
            // Shortcut method that is slightly more effecient, but combines with AdvanceNeuronState as well
            //State.Phase = 2 * Math.Atan(Alpha * spikeStats.Weight + c2);

            return new List<Spike>();
        }

        public override IEnumerable<Spike> AdvanceNeuronState(double newTime)
        {
            var outputSpikeTimes = new List<Spike>();
            var deltaTime = newTime - State.Time;
            var initialState = State;
            var c2 = GetC2(deltaTime);

            // Check for firing on trajectory; If Io>=0, we have to check the elapsed time between input spikes
            // and compare that to the unforced output spike frequency.
            if (Io > 0)
            {
                var thm = 2 * Math.Atan(c2);
                if (thm < initialState.StateVariable || deltaTime > Period)
                {
                    outputSpikeTimes.Add(new Spike {NeuronIndex = NeuronIndex , Time = NextSpikeTime});
                    // Check for case where phase completes multiple
                    // revolutions around phase circle between input spikes
                    if (deltaTime > Period)
                    {
                        var revCount = Math.Floor(deltaTime/Period);
                        for (var k = 2; k < revCount; k++)
                        {
                            outputSpikeTimes.Add(new Spike { NeuronIndex = NeuronIndex, Time = outputSpikeTimes.Last().Time + Period });
                        }
                        if (outputSpikeTimes.Last().Time + Period < newTime)
                        {
                            outputSpikeTimes.Add(new Spike { NeuronIndex = NeuronIndex, Time = outputSpikeTimes.Last().Time + Period });
                        }
                    }
                }
                State.StateVariable = thm;
            }
                // Check for firing on trajectory; If Io<0 the previous thp must be
                // between the positive fixed point and pi and thm must be between 
                // -pi and the negative fixed point in order for the neuron to have 
                // fired on trajectory, regardless of the weight.
            else
            {
                if (RemainingTimeToSpike < deltaTime)
                {
                    outputSpikeTimes.Add(new Spike { NeuronIndex = NeuronIndex, Time = NextSpikeTime });
                }                
                State.StateVariable = 2*Math.Atan(c2);
                
                if (outputSpikeTimes.Count > 0)
                {
                    State.StateVariable = PostSpikeState;
                }
            }

            State.Time = newTime;
            return outputSpikeTimes;
        }

        private double GetC2(double deltaTime)
        {
            if (Io < 0)
            {
                return Beta * (C1 + Math.Tanh(-Beta * deltaTime)) / (1 + C1 * Math.Tanh(-Beta * deltaTime));
            }
            return Beta * (C1 + Math.Tan(Beta * deltaTime)) / (1 - C1 * Math.Tan(Beta * deltaTime));
        }
    }
}