File size: 4,440 Bytes
cc38b25
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
/** Connectome LIF reference. Units: mV, ms, Hz. No fitted weights. */
export const PARAMETERS = Object.freeze({
  dt: 0.1,
  rest: -52,
  threshold: -45,
  tauM: 20,
  tauS: 5,
  refractory: 22,
  delay: 18,
  synapse: 0.275,
  poissonWeight: 68.75,
});
export function randomWord(i, t, seed = 1) {
  let x = (Math.imul(i + 1, 747796405) ^ Math.imul(t + 1, 2891336453) ^ seed) >>> 0;
  x = Math.imul(x ^ (x >>> 16), 2246822519) >>> 0;
  x = Math.imul(x ^ (x >>> 13), 3266489917) >>> 0;
  return (x ^ (x >>> 16)) >>> 0;
}
export function outgoingGraph(g) {
  const n = g.n,
    offsets = new Uint32Array(n + 1);
  for (const i of g.sources) offsets[i + 1]++;
  for (let i = 0; i < n; i++) offsets[i + 1] += offsets[i];
  const cursor = offsets.slice(),
    targets = new Uint32Array(g.sources.length),
    counts = new Uint32Array(g.sources.length);
  for (let j = 0; j < n; j++)
    for (let e = g.offsets[j]; e < g.offsets[j + 1]; e++) {
      const k = cursor[g.sources[e]]++;
      targets[k] = j;
      counts[k] = g.counts[e];
    }
  return { offsets, targets, counts };
}
const EM = Math.exp(-PARAMETERS.dt / PARAMETERS.tauM);
const ES = Math.exp(-PARAMETERS.dt / PARAMETERS.tauS);
const COUPLING = (PARAMETERS.tauS / (PARAMETERS.tauM - PARAMETERS.tauS)) * (EM - ES);

/** Event-driven scheduling; exactly resting neurons need no state update. */
export class BrainCPU {
  constructor(graph, { seed = 1 } = {}) {
    this.graph = graph;
    this.n = graph.n;
    this.seed = seed;
    this.out = outgoingGraph(graph);
    this.reset();
  }
  reset() {
    const n = this.n;
    this.v = new Float32Array(n).fill(-52);
    this.g = new Float32Array(n);
    this.until = new Uint32Array(n);
    this.counts = new Uint32Array(n);
    this.history = Array.from({ length: 19 }, () => []);
    this.tick = 0;
    this.active = new Uint32Array(n);
    this.present = new Uint8Array(n);
    this.activeCount = 0;
  }
  activate(i) {
    if (!this.present[i]) {
      this.present[i] = 1;
      this.active[this.activeCount++] = i;
    }
  }
  step(rates, externalEvents = null, silenced = false) {
    const driven = [];
    for (let i = 0; i < this.n; i++)
      if (externalEvents ? externalEvents[i] : rates[i] > 0) driven.push(i);
    return this.advance(rates, driven, externalEvents, silenced);
  }
  advance(rates, driven, externalEvents, silenced) {
    const { v, g, until, counts, active, present } = this,
      t = this.tick,
      p = PARAMETERS;
    const fired = [];
    let kept = 0;
    for (let k = 0; k < this.activeCount; k++) {
      const i = active[k];
      // No epsilon cutoff: skip only the exact stationary state. Refractory
      // deadlines remain stored and incoming events still check them.
      if (v[i] === p.rest && g[i] === 0) {
        present[i] = 0;
        continue;
      }
      active[kept++] = i;
      if (t >= until[i]) {
        v[i] = p.rest + (v[i] - p.rest) * EM + g[i] * COUPLING;
        g[i] *= ES;
        if (v[i] > p.threshold) fired.push(i);
      }
    }
    this.activeCount = kept;
    // Preserve the dense reference's source order and Float32 accumulation.
    fired.sort((a, b) => a - b);
    const due = this.history[t % 19],
      out = this.out;
    if (!silenced)
      for (const i of due) {
        const sign = this.graph.sign[i] * p.synapse;
        if (!sign) continue;
        for (let e = out.offsets[i]; e < out.offsets[i + 1]; e++) {
          const j = out.targets[e];
          if (t >= until[j]) {
            g[j] += out.counts[e] * sign;
            this.activate(j);
          }
        }
      }
    for (const i of driven) {
      if (
        t >= until[i] &&
        (externalEvents || randomWord(i, t, this.seed) / 4294967296 < (rates[i] * p.dt) / 1000)
      ) {
        v[i] += p.poissonWeight;
        this.activate(i);
      }
    }
    for (const i of fired) {
      v[i] = p.rest;
      g[i] = 0;
      until[i] = t + (rates[i] > 0 ? 0 : p.refractory);
      counts[i]++;
    }
    this.history[t % 19] = [];
    this.history[(t + p.delay) % 19] = fired;
    this.tick++;
    return fired;
  }
  batch(steps, rates, silenced = false) {
    const driven = [];
    for (let i = 0; i < this.n; i++) if (rates[i] > 0) driven.push(i);
    this.counts.fill(0);
    let total = 0;
    for (let k = 0; k < steps; k++) total += this.advance(rates, driven, null, silenced).length;
    return { counts: this.counts.slice(), tick: this.tick, total };
  }
}