change to simulating in potential space, rather than phase space

src/base.js

@@ -226,9 +226,9 @@ export("Network", Network);
 
 
 function Neuron(options) {
-  this.initial_phase = options.initial_phase != undefined ? options.initial_phase : Math.random();
   this.gamma = options.gamma;
   this.I = options.I;
+  this.initial_phase = this.g(options.initial_phase != undefined ? options.initial_phase : Math.random());
   this.connections = [];
 };
 
@@ -237,21 +237,29 @@ Neuron.prototype.initialize = function(current_time) {
 }
 
 Neuron.prototype.reset = function(current_time) {
-  this.set_phase(current_time, 0);
+  this.set_potential(current_time, 0);
 }
 
-Neuron.prototype.next_reset = function(current_time) {
-  return current_time + (1 - this.current_phase(current_time))*this.T();
+Neuron.prototype.next_reset = function() {
+  return this.last_potential.time + Math.log( (this.I-this.gamma*this.last_potential.potential)/(this.I-this.gamma)  );
 }
 
 Neuron.prototype.current_phase = function(current_time) {
-  return (current_time - this.last_spike.time)/this.T() + this.last_spike.phase;
+  return this.g(this.current_potential(current_time));
 }
 
-Neuron.prototype.set_phase = function(current_time, phase) {
-  this.last_spike = {time: current_time, phase: Math.max(0, phase)};
+Neuron.prototype.current_potential = function(current_time) {
+  return (1/this.gamma)*(this.I-(this.I-this.gamma*this.last_potential.potential)*Math.exp(-this.gamma*(current_time-this.last_potential.time)) );
+}
+
+Neuron.prototype.set_potential = function(current_time, potential) {
+  this.last_potential = {time: current_time, potential: potential};
 }
 
+Neuron.prototype.set_phase = function(current_time, phase) {
+  this.set_potential(current_time, this.f(phase));
+};
+
 Neuron.prototype.connect = function(post_synaptic, delay, strength, label) {
   this.connections.push({
     pre_synaptic: this,
@@ -262,23 +270,22 @@ Neuron.prototype.connect = function(post_synaptic, delay, strength, label) {
   });
 }
 
-Neuron.prototype.update_phase = function(current_time, new_phase) {
-  if(isNaN(new_phase) || new_phase > 1) {
-    this.set_phase(current_time, 0);
+Neuron.prototype.update_potential = function(current_time, new_potential) {
+  if(isNaN(new_potential) || new_potential > 1) {
+    this.set_potential(current_time, 0);
     return true;
   } else {
-    this.set_phase(current_time, Math.max(0, new_phase));
+    this.set_potential(current_time, new_potential);
     return false;
   }
 }
 
 Neuron.prototype.receive_spike = function(current_time, strength) {
-  var new_phase = this.phase_jump(this.current_phase(current_time), strength);
-  return this.update_phase(current_time, new_phase);
+  return this.update_potential(current_time, strength+this.current_potential(current_time));
 }
 
 Neuron.prototype.receive_phase_shift = function(current_time, phase_shift) {
-  return this.update_phase(current_time, this.current_phase(current_time) + phase_shift);
+  return this.update_potential(current_time, this.f(this.current_phase(current_time)+phase_shift));
 }
 
 Neuron.prototype.phase_jump = function(current_phase, strength) {

tests/test_base.js

@@ -84,7 +84,7 @@ var standard_period = 3.2580965380214812;
 test("initializing", function() {
   neuron.initial_phase = 0.5;
   neuron.initialize(0);
-  equals(neuron.current_phase(0), 0.5);
+  almost_equals(neuron.current_phase(0), 0.5);
 });
 
 test("calculating time scale", function() {
@@ -100,7 +100,7 @@ test("receiving reset sets phase to zero", function() {
 
 test("receiving spike will cause phase jump", function() {
   neuron.receive_spike(0.2*neuron.T(), 0.1);
-  equals(neuron.current_phase(0.2*neuron.T()), 0.26259348106668534);
+  almost_equals(neuron.current_phase(0.2*neuron.T()), 0.26259348106668534);
 });
 
 test("receiving reset after spike will reset phase to zero", function() {
@@ -109,18 +109,6 @@ test("receiving reset after spike will reset phase to zero", function() {
   almost_equals(neuron.current_phase(0.8+0.2*neuron.T()), 0.2);
 });
 
-test("receiving spike will not set phase beyond 1", function() {
-  neuron.receive_spike(0.9*neuron.T(), 0.5);
-  almost_equals(neuron.current_phase(0.9*neuron.T()), 0);
-  almost_equals(neuron.current_phase(neuron.T()), 0.1);
-});
-
-test("receiving spike will not set phase below 0", function() {
-  neuron.receive_spike(0.1*neuron.T(), -0.5);
-  almost_equals(neuron.current_phase(0.1*neuron.T()), 0);
-  almost_equals(neuron.current_phase(0.2*neuron.T()), 0.1);
-});
-
 test("receiving sub threshold spike will not signal reset", function() {
   ok(!neuron.receive_spike(0.1, 0.1));
 });
@@ -129,6 +117,13 @@ test("receiving supra threshold spike will signal reset", function() {
   ok(neuron.receive_spike(0.9, 0.5));
 });
 
+test("calculating time of next reset", function() {
+  almost_equals(neuron.next_reset(0), neuron.T());
+  almost_equals(neuron.next_reset(neuron.T()/2), neuron.T());
+  neuron.set_phase(0, 0.5);
+  almost_equals(neuron.next_reset(0), neuron.T()/2);
+});
+
 test("connecting to other neurons", function() {
   var n2 = new Neuron({ C: 1.04, gamma: 1 });
   neuron.connect(n2, 0.3, 0.4);
@@ -446,7 +441,7 @@ asyncTest("simulating dynamics of single transmitted spike", function() {
   simulator.start(0.5*n1.T(), function() {
     equals(n1.current_phase(0.5*n1.T()), 0.4);
     almost_equals(n2.current_phase(0.5*n1.T()), 0.5608294484932627);
-    almost_equals(n2.last_spike.time, 0.1*n1.T()+0.3);
+    almost_equals(n2.last_potential.time, 0.1*n1.T()+0.3);
     equals(simulator.past_events.length, 3);
 
     start();
@@ -478,7 +473,7 @@ asyncTest("simulating dynamics of two consequtive spikes sent to one neuron", fu
   simulator.start(1.5*n1.T(), function() {
     almost_equals(network.neurons[0].current_phase(1.5*n1.T()), 0.4);
     almost_equals(network.neurons[1].current_phase(1.5*n1.T()), 0.8897651188987299);
-    almost_equals(network.neurons[1].last_spike.time, 1.1*n1.T()+0.3);
+    almost_equals(network.neurons[1].last_potential.time, 1.1*n1.T()+0.3);
     equals(simulator.past_events.length, 6);
 
     start();
@@ -498,7 +493,7 @@ asyncTest("simulating dynamics of two neurons in different sub networks", functi
   simulator.start(1.5*n1.T(), function() {
     almost_equals(net1.neurons[0].current_phase(1.5*n1.T()), 0.4);
     almost_equals(net2.neurons[0].current_phase(1.5*n1.T()), 0.8897651188987299);
-    almost_equals(net2.neurons[0].last_spike.time, 1.1*n1.T()+0.3);
+    almost_equals(net2.neurons[0].last_potential.time, 1.1*n1.T()+0.3);
     equals(simulator.past_events.length, 6);
 
     start();