Coverage for src/driada/experiment/synthetic/manifold_circular.py: 94.52%
73 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-07-25 15:40 +0300
1"""
2Circular manifold generation for head direction cells.
4This module contains functions for generating synthetic neural data on circular
5manifolds, typically used to model head direction cells.
6"""
8import numpy as np
9from .core import validate_peak_rate, generate_pseudo_calcium_signal
10from .utils import get_effective_decay_time
11from ..exp_base import Experiment
12from ...information.info_base import TimeSeries, MultiTimeSeries
15def generate_circular_random_walk(length, step_std=0.1, seed=None):
16 """
17 Generate a random walk on a circle (head direction trajectory).
19 Parameters
20 ----------
21 length : int
22 Number of time points.
23 step_std : float
24 Standard deviation of angular steps in radians.
25 seed : int, optional
26 Random seed for reproducibility.
28 Returns
29 -------
30 angles : ndarray
31 Array of angles in radians [0, 2π).
32 """
33 if seed is not None:
34 np.random.seed(seed)
36 # Generate random steps
37 steps = np.random.normal(0, step_std, length)
39 # Cumulative sum to get trajectory
40 angles = np.cumsum(steps)
42 # Wrap to [0, 2π)
43 angles = angles % (2 * np.pi)
45 return angles
48def von_mises_tuning_curve(angles, preferred_direction, kappa):
49 """
50 Calculate neural response using Von Mises tuning curve.
52 Parameters
53 ----------
54 angles : ndarray
55 Current head directions in radians.
56 preferred_direction : float
57 Preferred direction of the neuron in radians.
58 kappa : float
59 Concentration parameter (inverse width of tuning curve).
60 Higher kappa = narrower tuning.
62 Returns
63 -------
64 response : ndarray
65 Neural response (firing rate modulation).
66 """
67 # Von Mises distribution normalized to max=1
68 response = np.exp(kappa * (np.cos(angles - preferred_direction) - 1))
69 return response
72def generate_circular_manifold_neurons(n_neurons, head_direction, kappa=4.0,
73 baseline_rate=0.1, peak_rate=1.0,
74 noise_std=0.05, seed=None):
75 """
76 Generate population of head direction cells with Von Mises tuning.
78 Parameters
79 ----------
80 n_neurons : int
81 Number of neurons in the population.
82 head_direction : ndarray
83 Head direction trajectory in radians.
84 kappa : float
85 Concentration parameter for Von Mises tuning curves.
86 Typical values: 2-8 (higher = narrower tuning).
87 baseline_rate : float
88 Baseline firing rate when far from preferred direction.
89 Default is 0.1 Hz (realistic for sparse firing neurons).
90 peak_rate : float
91 Peak firing rate at preferred direction.
92 Default is 1.0 Hz (realistic for calcium imaging).
93 Values >2 Hz may cause calcium signal saturation.
94 noise_std : float
95 Standard deviation of noise in firing rates.
96 seed : int, optional
97 Random seed for reproducibility.
99 Returns
100 -------
101 firing_rates : ndarray
102 Shape (n_neurons, n_timepoints) with firing rates.
103 preferred_directions : ndarray
104 Preferred direction for each neuron in radians.
105 """
106 # Validate firing rate
107 validate_peak_rate(peak_rate, context="generate_circular_manifold_neurons")
109 if seed is not None:
110 np.random.seed(seed)
112 n_timepoints = len(head_direction)
114 # Uniformly distribute preferred directions around the circle
115 preferred_directions = np.linspace(0, 2*np.pi, n_neurons, endpoint=False)
117 # Add small random jitter to break perfect symmetry
118 jitter = np.random.normal(0, 0.1, n_neurons)
119 preferred_directions = (preferred_directions + jitter) % (2*np.pi)
121 # Generate firing rates for each neuron
122 firing_rates = np.zeros((n_neurons, n_timepoints))
124 for i in range(n_neurons):
125 # Von Mises tuning curve
126 tuning_response = von_mises_tuning_curve(head_direction,
127 preferred_directions[i],
128 kappa)
130 # Scale to desired firing rate range
131 firing_rate = baseline_rate + (peak_rate - baseline_rate) * tuning_response
133 # Add noise
134 noise = np.random.normal(0, noise_std, n_timepoints)
135 firing_rate = np.maximum(0, firing_rate + noise) # Ensure non-negative
137 firing_rates[i, :] = firing_rate
139 return firing_rates, preferred_directions
142def generate_circular_manifold_data(n_neurons, duration=600, sampling_rate=20.0,
143 kappa=4.0, step_std=0.1,
144 baseline_rate=0.1, peak_rate=1.0,
145 noise_std=0.05,
146 decay_time=2.0, calcium_noise_std=0.1,
147 seed=None, verbose=True):
148 """
149 Generate synthetic data with neurons on circular manifold (head direction cells).
151 Parameters
152 ----------
153 n_neurons : int
154 Number of neurons.
155 duration : float
156 Duration in seconds.
157 sampling_rate : float
158 Sampling rate in Hz.
159 kappa : float
160 Von Mises concentration parameter (tuning width).
161 step_std : float
162 Standard deviation of head direction random walk steps.
163 baseline_rate : float
164 Baseline firing rate. Default is 0.1 Hz.
165 peak_rate : float
166 Peak firing rate at preferred direction. Default is 1.0 Hz.
167 Values >2 Hz may cause calcium signal saturation.
168 noise_std : float
169 Noise in firing rates.
170 decay_time : float
171 Calcium decay time constant.
172 calcium_noise_std : float
173 Noise in calcium signal.
174 seed : int, optional
175 Random seed.
176 verbose : bool
177 Print progress.
179 Returns
180 -------
181 calcium_signals : ndarray
182 Calcium signals (n_neurons x n_timepoints).
183 head_direction : ndarray
184 Head direction trajectory.
185 preferred_directions : ndarray
186 Preferred direction for each neuron.
187 firing_rates : ndarray
188 Underlying firing rates.
189 """
190 if seed is not None:
191 np.random.seed(seed)
193 n_timepoints = int(duration * sampling_rate)
195 if verbose:
196 print(f'Generating circular manifold data: {n_neurons} neurons, {duration}s')
198 # Generate head direction trajectory
199 if verbose:
200 print(' Generating head direction trajectory...')
201 head_direction = generate_circular_random_walk(n_timepoints, step_std, seed)
203 # Generate neural responses
204 if verbose:
205 print(' Generating neural responses with Von Mises tuning...')
206 firing_rates, preferred_directions = generate_circular_manifold_neurons(
207 n_neurons, head_direction, kappa,
208 baseline_rate, peak_rate, noise_std,
209 seed=(seed + 1) if seed else None
210 )
212 # Convert firing rates to calcium signals
213 if verbose:
214 print(' Converting to calcium signals...')
215 calcium_signals = np.zeros((n_neurons, n_timepoints))
217 for i in range(n_neurons):
218 # Generate Poisson events from firing rates
219 prob_spike = firing_rates[i, :] / sampling_rate
220 prob_spike = np.clip(prob_spike, 0, 1) # Ensure valid probability
221 events = np.random.binomial(1, prob_spike)
223 # Convert to calcium using existing function
224 calcium_signal = generate_pseudo_calcium_signal(
225 events=events,
226 duration=duration,
227 sampling_rate=sampling_rate,
228 amplitude_range=(0.5, 2.0),
229 decay_time=decay_time,
230 noise_std=calcium_noise_std
231 )
232 calcium_signals[i, :] = calcium_signal
234 if verbose:
235 print(' Done!')
237 return calcium_signals, head_direction, preferred_directions, firing_rates
240def generate_circular_manifold_exp(n_neurons=100, duration=600, fps=20.0,
241 kappa=4.0, step_std=0.1,
242 baseline_rate=0.1, peak_rate=1.0,
243 noise_std=0.05,
244 decay_time=2.0, calcium_noise_std=0.1,
245 add_mixed_features=False,
246 seed=None, verbose=True, return_info=False):
247 """
248 Generate complete experiment with circular manifold (head direction cells).
250 Parameters
251 ----------
252 n_neurons : int
253 Number of neurons.
254 duration : float
255 Duration in seconds.
256 fps : float
257 Sampling rate (frames per second).
258 kappa : float
259 Von Mises concentration parameter.
260 step_std : float
261 Head direction random walk step size.
262 baseline_rate : float
263 Baseline firing rate. Default is 0.1 Hz.
264 peak_rate : float
265 Peak firing rate. Default is 1.0 Hz.
266 noise_std : float
267 Firing rate noise.
268 decay_time : float
269 Calcium decay time.
270 calcium_noise_std : float
271 Calcium signal noise.
272 add_mixed_features : bool
273 Whether to add circular_angle MultiTimeSeries (cos/sin representation).
274 seed : int, optional
275 Random seed.
276 verbose : bool
277 Print progress.
279 Returns
280 -------
281 exp : Experiment
282 DRIADA Experiment object with circular manifold data.
283 """
284 # Calculate effective decay time for shuffle mask
285 effective_decay_time = get_effective_decay_time(decay_time, duration, verbose)
287 # Generate data
288 calcium, head_direction, preferred_directions, firing_rates = generate_circular_manifold_data(
289 n_neurons=n_neurons,
290 duration=duration,
291 sampling_rate=fps,
292 kappa=kappa,
293 step_std=step_std,
294 baseline_rate=baseline_rate,
295 peak_rate=peak_rate,
296 noise_std=noise_std,
297 decay_time=decay_time,
298 calcium_noise_std=calcium_noise_std,
299 seed=seed,
300 verbose=verbose
301 )
303 # Create static features
304 static_features = {
305 'fps': fps,
306 't_rise_sec': 0.04,
307 't_off_sec': effective_decay_time, # Use effective decay time for shuffle mask
308 'manifold_type': 'circular',
309 'kappa': kappa,
310 'baseline_rate': baseline_rate,
311 'peak_rate': peak_rate,
312 }
314 # Create dynamic features
315 head_direction_ts = TimeSeries(
316 data=head_direction,
317 discrete=False
318 )
320 dynamic_features = {
321 'head_direction': head_direction_ts
322 }
324 # Add circular_angle MultiTimeSeries if requested
325 if add_mixed_features:
326 # Create circular_angle as MultiTimeSeries with cos and sin components
327 # This is the proper representation for circular variables
328 cos_component = np.cos(head_direction)
329 sin_component = np.sin(head_direction)
330 circular_angle_mts = MultiTimeSeries([
331 TimeSeries(data=cos_component, discrete=False),
332 TimeSeries(data=sin_component, discrete=False)
333 ])
334 dynamic_features['circular_angle'] = circular_angle_mts
336 # Store additional information
337 static_features['preferred_directions'] = preferred_directions
339 # Create experiment
340 exp = Experiment(
341 signature='circular_manifold_exp',
342 calcium=calcium,
343 spikes=None, # Will be extracted from calcium if needed
344 static_features=static_features,
345 dynamic_features=dynamic_features,
346 exp_identificators={
347 'manifold': 'circular',
348 'n_neurons': n_neurons,
349 'duration': duration
350 }
351 )
353 # Store firing rates as additional data
354 exp.firing_rates = firing_rates
356 # Create info dictionary if requested
357 if return_info:
358 info = {
359 'manifold_type': 'circular',
360 'n_neurons': n_neurons,
361 'head_direction': head_direction,
362 'preferred_directions': preferred_directions,
363 'firing_rates': firing_rates,
364 'parameters': {
365 'kappa': kappa,
366 'step_std': step_std,
367 'baseline_rate': baseline_rate,
368 'peak_rate': peak_rate,
369 'noise_std': noise_std,
370 'decay_time': decay_time,
371 'calcium_noise_std': calcium_noise_std
372 }
373 }
374 return exp, info
376 return exp