Add Kuramoto example

nickmcintyre · Oct 26, 2019 · 88688ef · 88688ef
1 parent 720ec15
commit 88688ef
Show file tree

Hide file tree

Showing 8 changed files with 439 additions and 12 deletions.
diff --git a/README.md b/README.md
@@ -1,4 +1,5 @@
 # número
+![número](numero.png)
 > A friendly and intuitive math library for p5.js
 
 This project aims to bring the mathematical chops of [NumPy](https://numpy.org/) and [SymPy](https://www.sympy.org/en/index.html) to the [p5.js](https://p5js.org/) ecosystem. A few guiding principles are:
@@ -16,5 +17,17 @@ Excellent libraries like [math.js](https://mathjs.org/) and [p5.dimensions](http
 
 Also, math + code = awesome :)
 
+## Usage
+
+```javascript
+const b = num.tidy(() => {
+    const a = createTensor([[1, 2], [3, 4]]);
+    const x = createTensor([5, 6]);
+    return a.dot(x);
+});
+
+print(b.toString());
+```
+
 ## Contributing
 See [CONTRIBUTING](CONTRIBUTING.md).
diff --git a/dist/numero.js b/dist/numero.js
diff --git a/examples/kuramoto/index.html b/examples/kuramoto/index.html
@@ -0,0 +1,10 @@
+<html>
+  <head>
+    <title>Kuramoto Example</title>
+    <script src="https://cdnjs.cloudflare.com/ajax/libs/p5.js/0.9.0/p5.js"></script>
+    <script src="../../dist/numero.js"></script>
+    <script src="sketch.js"></script>
+  </head>
+  <body>
+  </body>
+</html>
diff --git a/examples/kuramoto/license.txt b/examples/kuramoto/license.txt
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2019 Nick McIntyre
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/examples/kuramoto/sketch.js b/examples/kuramoto/sketch.js
@@ -0,0 +1,259 @@
+// Simulate the synchronization of neurons in the brain by varying their connective
+// arrangement.
+// https://researchspace.auckland.ac.nz/bitstream/handle/2292/2666/esc-tr-638-1.pdf
+const networkSize = 10;
+let time = 0;
+const dt = 0.01;
+let coupling;
+const couplingStrength = 5;
+const noiseLevel = 0.5;
+let naturalFrequency;
+let phase;
+let velocity;
+let acceleration;
+let arrangement = 'Press a key 1-5';
+
+function setup() {
+  createCanvas(400, 400);
+  naturalFrequency = num.Tensor.random([networkSize]).mult(TWO_PI);
+  phase = num.Tensor.zeros([networkSize]);
+  velocity = num.Tensor.zeros([networkSize]);
+  acceleration = num.Tensor.zeros([networkSize]);
+  coupling = num.Tensor.zeros([networkSize, networkSize]);
+}
+
+function draw() {
+  background(255);
+  noStroke();
+  fill(54, 86, 148);
+  text(arrangement, 10, 20);
+  translate(width / 2, height / 2);
+  const theta = phase.arraySync();
+  for (let i = 0; i < networkSize; i += 1) {
+    const r = i * TWO_PI / networkSize;
+    const a = map(theta[i], 0, TWO_PI, 0, 255);
+    push();
+    rotate(r);
+    translate(100, 0);
+    fill(54, 86, 148, a);
+    circle(0, 0, 25);
+    pop();
+  }
+
+  step();
+}
+
+function keyPressed() {
+  switch (key) {
+    case '1': {
+      /**
+       * Linear Unidirectional
+       *
+       * [[0, X, 0, 0],
+       *  [0, 0, X, 0],
+       *  [0, 0, 0, X],
+       *  [0, 0, 0, 0]]
+       */
+      coupling.dispose();
+      const cpl = new Array(networkSize);
+      for (let i = 0; i < networkSize; i += 1) {
+        cpl[i] = new Array(networkSize);
+        for (let j = 0; j < networkSize; j += 1) {
+          if (j === i + 1) {
+            cpl[i][j] = couplingStrength;
+          } else {
+            cpl[i][j] = 0;
+          }
+        }
+      }
+
+      coupling = createTensor(cpl);
+      arrangement = 'Linear Unidirectional';
+      break;
+    }
+    case '2': {
+      /**
+       * Linear Bidirectional
+       *
+       * [[0, X, 0, 0],
+       *  [X, 0, X, 0],
+       *  [0, X, 0, X],
+       *  [0, 0, X, 0]]
+       */
+      coupling.dispose();
+      const cpl = new Array(networkSize);
+      for (let i = 0; i < networkSize; i += 1) {
+        cpl[i] = new Array(networkSize);
+        for (let j = 0; j < networkSize; j += 1) {
+          if (j === i + 1) {
+            cpl[i][j] = couplingStrength;
+          } else if (i === j + 1) {
+            cpl[i][j] = couplingStrength;
+          } else {
+            cpl[i][j] = 0;
+          }
+        }
+      }
+
+      coupling = createTensor(cpl);
+      arrangement = 'Linear Bidirectional';
+      break;
+    }
+    case '3': {
+      /**
+       * Box Unidirectional
+       *
+       * [[0, X, 0, 0],
+       *  [0, 0, X, 0],
+       *  [0, 0, 0, X],
+       *  [X, 0, 0, 0]]
+       */
+      coupling.dispose();
+      const cpl = new Array(networkSize);
+      for (let i = 0; i < networkSize; i += 1) {
+        cpl[i] = new Array(networkSize);
+        for (let j = 0; j < networkSize; j += 1) {
+          if (j === i + 1) {
+            cpl[i][j] = couplingStrength;
+          } else if (j === 0 && i === networkSize - 1) {
+            cpl[i][j] = couplingStrength;
+          } else {
+            cpl[i][j] = 0;
+          }
+        }
+      }
+
+      coupling = createTensor(cpl);
+      arrangement = 'Box Unidirectional';
+      break;
+    }
+    case '4': {
+      /**
+       * Box Bidirectional
+       *
+       * [[0, X, 0, X],
+       *  [X, 0, X, 0],
+       *  [0, X, 0, X],
+       *  [X, 0, X, 0]]
+       */
+      coupling.dispose();
+      const cpl = new Array(networkSize);
+      for (let i = 0; i < networkSize; i += 1) {
+        cpl[i] = new Array(networkSize);
+        for (let j = 0; j < networkSize; j += 1) {
+          if (j === i + 1) {
+            cpl[i][j] = couplingStrength;
+          } else if (i === j + 1) {
+            cpl[i][j] = couplingStrength;
+          } else if (j === 0 && i === networkSize - 1) {
+            cpl[i][j] = couplingStrength;
+          } else if (i === 0 && j === networkSize - 1) {
+            cpl[i][j] = couplingStrength;
+          } else {
+            cpl[i][j] = 0;
+          }
+        }
+      }
+
+      coupling = createTensor(cpl);
+      arrangement = 'Box Bidirectional';
+      break;
+    }
+    case '5': {
+      /**
+       * All-to-all
+       *
+       * [[0, X, X, X],
+       *  [X, 0, X, X],
+       *  [X, X, 0, X],
+       *  [X, X, X, 0]]
+       */
+      coupling.dispose();
+      const cpl = new Array(networkSize);
+      for (let i = 0; i < networkSize; i += 1) {
+        cpl[i] = new Array(networkSize);
+        for (let j = 0; j < networkSize; j += 1) {
+          if (i === j) {
+            cpl[i][j] = 0;
+          } else {
+            cpl[i][j] = couplingStrength;
+          }
+        }
+      }
+
+      coupling = createTensor(cpl);
+      arrangement = 'All-to-All';
+      break;
+    }
+    default: {
+      break;
+    }
+  }
+}
+
+/**
+ * Calculate the next phase of the oscillator network by solving the governing
+ * equation with the classical Runge-Kutta method.
+ *
+ * https://en.wikipedia.org/wiki/Runge-Kutta_methods
+ */
+function step() {
+  time += dt;
+  const zeta = noiseLevel * noise(time);
+  const result = num.tidy(() => {
+    const oldPhase = phase.copy();
+    const oldVelocity = velocity.copy();
+
+    const zeros = num.Tensor.zeros([networkSize]);
+    const k1 = diff(zeros, zeta).mult(dt);
+    const k2 = diff(k1.div(2), zeta).mult(dt);
+    const k3 = diff(k2.div(2), zeta).mult(dt);
+    const k4 = diff(k3, zeta).mult(dt);
+    const solution = k1.add(k2.mult(2)).add(k3.mult(2)).add(k4).div(6);
+
+    const newPhase = phase.add(solution).mod(TWO_PI);
+    const newVelocity = phase.sub(oldPhase).div(dt);
+    const newAcceleration = velocity.sub(oldVelocity).div(dt);
+
+    return {
+      newPhase,
+      newVelocity,
+      newAcceleration,
+    };
+  });
+
+  phase.dispose();
+  velocity.dispose();
+  acceleration.dispose();
+
+  phase = result.newPhase;
+  velocity = result.newVelocity;
+  acceleration = result.newAcceleration;
+}
+
+/**
+ * Calculate the time derivative of an oscillator's phase using the Kuramoto model.
+ *
+ * https://en.wikipedia.org/wiki/Kuramoto_model
+ */
+function diff(increment, zeta) {
+  const t = num.tidy(() => {
+    let dTheta = new Array(networkSize);
+    for (let i = 0; i < networkSize; i += 1) {
+      dTheta[i] = phase.sub(increment);
+    }
+
+    dTheta = num.Tensor.stack(dTheta);
+    const t_ = dTheta.sin()
+                      .mult(coupling)
+                      .sum(0)
+                      .div(networkSize)
+                      .add(zeta)
+                      .add(naturalFrequency);
+
+    return t_;
+  });
+  const result = createTensor(t);
+
+  return result;
+}
diff --git a/numero.png b/numero.png