Skip to content

Commit bab424c

Browse files
sebsjamessebsjames
authored
Merge pull request #134 from sebsjames/dev/random_walk
A random walk class after Stone et al
2 parents e5c1bde + 7493784 commit bab424c

3 files changed

Lines changed: 237 additions & 0 deletions

File tree

‎cmake/module_definitions.cmake‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,10 @@ macro(setup_module_variables_for_maths base_directory json_directory)
6565
${SM_MAT_MODULES}
6666
${base_directory}/sm/spline.cppm
6767
)
68+
set(SM_RANDOM_WALK_MODULES
69+
${SM_SPLINE_MODULES}
70+
${base_directory}/sm/random_walk.cppm
71+
)
6872
set(SM_BEZCOORD_MODULES
6973
${SM_VEC_MODULES}
7074
${base_directory}/sm/bezcoord.cppm

‎sm/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ install(
3838
polysolve.cppm
3939
quaternion.cppm
4040
random.cppm
41+
random_walk.cppm
4142
range.cppm
4243
rect.cppm
4344
scale.cppm

‎sm/random_walk.cppm‎

Lines changed: 232 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,232 @@
1+
// -*- C++ -*-
2+
/*!
3+
* This file is part of sebsjames/maths, a library of maths code for modern C++
4+
*
5+
* See https://github.com/sebsjames/maths
6+
*
7+
* \file
8+
*
9+
* A Random walk module. Generate a bee-like random walk, following the description given in:
10+
*
11+
* "An Anatomically Constrained Model for Path Integration in the Bee
12+
* Brain", Current Biology 27, 3069-3085, October 23, 2017 (Stone et al.)
13+
*
14+
* DOI: http://dx.doi.org/10.1016/j.cub.2017.08.052
15+
*
16+
* \author Seb James
17+
* \date 2026
18+
*/
19+
module;
20+
21+
#include <cstdint>
22+
#include <memory>
23+
#include <cmath>
24+
#include <stdexcept>
25+
26+
export module sm.random_walk;
27+
28+
import sm.spline;
29+
import sm.mathconst;
30+
import sm.random;
31+
import sm.algo;
32+
import sm.vec;
33+
import sm.vvec;
34+
35+
export namespace sm
36+
{
37+
// Make a randomized path to follow
38+
template <typename T>
39+
struct random_walk
40+
{
41+
random_walk() {}
42+
43+
random_walk (const std::uint32_t _n_steps, const std::uint32_t _a_tau)
44+
{
45+
this->n_steps = _n_steps;
46+
this->a_tau = _a_tau;
47+
this->init();
48+
}
49+
50+
random_walk (const std::uint32_t _n_steps, const std::uint32_t _a_tau, const T& _kappa)
51+
{
52+
this->n_steps = _n_steps;
53+
this->a_tau = _a_tau;
54+
this->kappa = _kappa;
55+
this->init();
56+
}
57+
58+
random_walk (const std::uint32_t _n_steps, const std::uint32_t _a_tau, const T& _kappa, const T& acc_max)
59+
{
60+
this->n_steps = _n_steps;
61+
this->a_tau = _a_tau;
62+
this->kappa = _kappa;
63+
this->amm.max = acc_max;
64+
this->init();
65+
}
66+
67+
void init()
68+
{
69+
this->rVM = std::make_unique<sm::rand_vonmises<T>> (T{0}, kappa);
70+
this->a.resize (this->n_steps / this->a_tau);
71+
this->a.randomize();
72+
this->a *= (amm.span());
73+
this->a += amm.min;
74+
cubic_spline_expansion (this->a, this->a_tau);
75+
}
76+
77+
// Reset state
78+
void reset()
79+
{
80+
this->t = 0;
81+
this->theta = T{0};
82+
this->omega = T{0};
83+
this->velocity = {};
84+
this->speed = T{0};
85+
}
86+
87+
void about_turn() { this->theta += sm::mathconst<float>::pi; }
88+
89+
// Advance the route generation by one timestep
90+
void step()
91+
{
92+
// This is the model as stated in the paper and it should be equivalent to lfilter
93+
// function.
94+
T epsilon = this->rVM->get(); // Angular acceleration
95+
this->omega = this->lambda * this->omega + epsilon;
96+
this->theta += this->omega;
97+
98+
T accel = T{0};
99+
if (t < this->a.size()) { accel = this->a[this->t]; }
100+
101+
sm::vec<T, 2> thrust = { accel * std::sin (theta), accel * std::cos (theta) };
102+
this->velocity = (this->velocity + thrust) * one_minus_FD;
103+
this->speed = (this->speed + accel) * one_minus_FD;
104+
105+
++this->t;
106+
if (this->t > this->n_steps) { this->reset(); }
107+
}
108+
109+
// Number of steps total.
110+
std::uint32_t n_steps = 0;
111+
// Current time step
112+
std::uint32_t t = 0;
113+
114+
// State
115+
T theta = T{0}; // Heading/theta
116+
T omega = T{0}; // Angular velocity
117+
sm::vec<T, 2> velocity = {}; // Cartesian velocity (or can use speed):
118+
T speed = T{0}; // Linear speed
119+
120+
// Parameters
121+
const T lambda = T{0.4};
122+
T kappa = T{100}; // Von Mises concentration parameter
123+
sm::vvec<T> a = {}; // Acceleration values
124+
// Uniform RNG range outbound [0, 0.05]
125+
sm::range<T> amm = { T{0}, T{0.005} };
126+
// how often does the acceleration change?
127+
std::uint32_t a_tau = 50;
128+
129+
// FD is the drag coefficient
130+
static constexpr T FD = T{0.15};
131+
static constexpr T one_minus_FD = (T{1} - FD);
132+
133+
// Random number generation
134+
std::unique_ptr<sm::rand_vonmises<T>> rVM;
135+
};
136+
137+
template <typename T, std::uint32_t N>
138+
void cubic_spline_expansion (sm::vvec<T>& v, std::uint32_t n)
139+
{
140+
if (v.size() != N) { throw std::runtime_error ("v.size != N!"); }
141+
sm::vec<sm::vec<T, 2>, N> p;
142+
for (std::uint32_t i = 0; i < N; ++i) { p[i] = { static_cast<T>(n) * i, v[i] }; }
143+
sm::spline<T, N> spl (p);
144+
v.resize (v.size() * (n + 1), T{0});
145+
T ti = T{0};
146+
for (std::uint32_t i = 0; i < v.size(); ++i, ti += T{1}) { v[i] = spl.compute (ti); }
147+
}
148+
149+
// Place n elements between each element in v, computing a cubic spline interpolation, assuming
150+
// that the x in f(x) increases by the value 1 with each step. Think of x as t, and there is 1
151+
// timestep between each element in the expanded vvec v. This calls fixed-size versions of
152+
// cubic_spline_expansion due to the limitation from our fixed size sm::mat.
153+
template <typename T>
154+
void cubic_spline_expansion (sm::vvec<T>& v, std::uint32_t n)
155+
{
156+
// sm::mat is a fixed size matrix template class, so we have to have a switch statement
157+
switch (v.size()) {
158+
case 2:
159+
cubic_spline_expansion<T, 2> (v, n);
160+
break;
161+
case 3:
162+
cubic_spline_expansion<T, 3> (v, n);
163+
break;
164+
case 4:
165+
cubic_spline_expansion<T, 4> (v, n);
166+
break;
167+
case 5:
168+
cubic_spline_expansion<T, 5> (v, n);
169+
break;
170+
case 6:
171+
cubic_spline_expansion<T, 6> (v, n);
172+
break;
173+
case 7:
174+
cubic_spline_expansion<T, 7> (v, n);
175+
break;
176+
case 8:
177+
cubic_spline_expansion<T, 8> (v, n);
178+
break;
179+
case 9:
180+
cubic_spline_expansion<T, 9> (v, n);
181+
break;
182+
case 10:
183+
cubic_spline_expansion<T, 10> (v, n);
184+
break;
185+
case 12:
186+
cubic_spline_expansion<T, 12> (v, n);
187+
break;
188+
case 14:
189+
cubic_spline_expansion<T, 14> (v, n);
190+
break;
191+
case 16:
192+
cubic_spline_expansion<T, 16> (v, n);
193+
break;
194+
case 18:
195+
cubic_spline_expansion<T, 18> (v, n);
196+
break;
197+
case 20:
198+
cubic_spline_expansion<T, 20> (v, n);
199+
break;
200+
case 30:
201+
cubic_spline_expansion<T, 30> (v, n);
202+
break;
203+
case 40:
204+
cubic_spline_expansion<T, 40> (v, n);
205+
break;
206+
case 50:
207+
cubic_spline_expansion<T, 50> (v, n);
208+
break;
209+
case 60:
210+
cubic_spline_expansion<T, 60> (v, n);
211+
break;
212+
case 70:
213+
cubic_spline_expansion<T, 70> (v, n);
214+
break;
215+
case 80:
216+
cubic_spline_expansion<T, 80> (v, n);
217+
break;
218+
case 90:
219+
cubic_spline_expansion<T, 90> (v, n);
220+
break;
221+
case 100:
222+
cubic_spline_expansion<T, 100> (v, n);
223+
break;
224+
case 0:
225+
case 1:
226+
default:
227+
throw std::runtime_error ("sm::cubic_spline_expansion cannot handle that sized input");
228+
break;
229+
}
230+
}
231+
232+
} // namespace

0 commit comments

Comments
 (0)