#!/usr/bin/env python3
# Filename: reservoir.py
# Author: Andrew Carr (direct questions to andrewcarr319@gmail.com after June 2026)
"""
A module to simulate the dynamics of a reservoir computer with an RK4
DDE solver along with a analysis class for visualizing and testing reservoirs.
It is built with JAX, specifically JIT compilation, for speed.
"""
OMP_NUM_THREADS=1
MKL_NUM_THREADS=1
OPENBLAS_NUM_THREADS=1
NUMEXPR_NUM_THREADS=1
import jax
import jax.numpy as jnp
from functools import partial
import numpy as np
import matplotlib.pyplot as plt
import os
from PIL import Image
from io import BytesIO
[docs]
class Reservoir:
def __init__(self, *, h=.02, theta, N, sd=0, amp=1, tau=0, fb_gain=0,
norm_factor=None, norm_offset=None, state_shape=1,
input_connectivity=0.2, mask=None, mask_seed=0,
load=False, normalize_mask=True, VDC=20, n_envs=1,
n_steps_needed=None):
self.h = h
self.theta = theta
self.N = N
self.sd = int(sd)
self.amplification = amp
self.tau = tau
self.fb_gain = fb_gain
self.mask_seed = mask_seed
self.VDC = VDC
self.tau_steps = int(self.tau * jnp.floor(self.theta/self.h))
m = self.tau_steps
self.buf_len = max(m+2, 2)
# init reservoir state even if will be reinited in algo later
# handles dqn and ppo
self._init_reservoir_state(n_envs)
self.stride = int(jnp.maximum(1, jnp.floor(self.theta / self.h)))
if self.sd > self.stride:
raise ValueError("Sample delay must be less than or equal to stride")
if n_steps_needed is None:
self.n_steps_needed = (self.N) * self.stride + 1
else:
self.n_steps_needed = n_steps_needed
# print(f"n_steps_needed: {self.n_steps_needed}")
self.normalization_factor = norm_factor
self.normalization_offset = norm_offset
if type(state_shape) is not int:
state_shape = int(state_shape[-1])
print(f"Note: state_shape was not int, using {state_shape}")
self.state_shape = state_shape
# print(state_shape)
self.input_connectivity = input_connectivity
if load and mask is not None:
self.mask = mask
else:
self.mask = self.create_mask(mask_seed, normalize_mask)
def _make_reservoir_state(self, n_envs: int):
return {
"pos_bufs": jnp.zeros((n_envs, self.buf_len), dtype=jnp.float32),
"vel_bufs": jnp.zeros((n_envs, self.buf_len), dtype=jnp.float32),
"buf_idxs": jnp.zeros((n_envs,), dtype=jnp.int32),
"buf_cnts": jnp.zeros((n_envs,), dtype=jnp.int32),
}
def _get_state(self, mode: str, env_idx: int):
state = self.res_state[mode]
return (
state["pos_bufs"][env_idx],
state["vel_bufs"][env_idx],
state["buf_idxs"][env_idx],
state["buf_cnts"][env_idx],
)
def _set_state(self, mode: str, env_idx: int, pos_buf, vel_buf, buf_idx, buf_cnt):
state = self.res_state[mode]
state["pos_bufs"] = state["pos_bufs"].at[env_idx].set(pos_buf)
state["vel_bufs"] = state["vel_bufs"].at[env_idx].set(vel_buf)
state["buf_idxs"] = state["buf_idxs"].at[env_idx].set(buf_idx)
state["buf_cnts"] = state["buf_cnts"].at[env_idx].set(buf_cnt)
self.res_state[mode] = state
def _init_reservoir_state(self, train_envs, val_envs=1, test_envs=1):
self.res_state = {
"train": self._make_reservoir_state(train_envs),
"val": self._make_reservoir_state(val_envs),
"test": self._make_reservoir_state(test_envs),
}
[docs]
def create_mask(self, seed, normalize_mask=False):
mask_rand = np.random.RandomState(seed)
mask = mask_rand.choice([-1,1],size=(self.state_shape, self.N))
# add sparsity to mask
REPLACE_COUNT = np.int32(self.input_connectivity * self.state_shape * self.N)
mask_reshaped = np.reshape(mask, [self.state_shape * self.N,])
mask_reshaped[mask_rand.choice(self.state_shape * self.N, REPLACE_COUNT, replace=False)] = 0
mask = np.reshape(mask_reshaped,[self.state_shape, self.N])
if normalize_mask:
mask = mask / self.state_shape
return mask
[docs]
def sim(self, *, obs=None, Vin=None, mode='train',env_idx=0,
full_data=False, direct_fb=False, get_Vin=False, update_state=True):
if Vin is None:
state = np.array(obs).reshape(1,self.state_shape)
u = np.divide(state + self.normalization_offset, self.normalization_factor)
Vin = np.dot(u, self.mask)
Vin = Vin.reshape(-1)
if get_Vin:
return Vin*self.amplification + self.VDC
# print(f'Vin mean new: {np.mean(Vin)}')
pos_buf, vel_buf, buf_idx, buf_cnt = self._get_state(mode, env_idx)
positions, theta_vals, use_fb, pos_buf, vel_buf, buf_idx, buf_cnt = Reservoir.rk4(
n_steps=self.n_steps_needed,
# y0=jnp.asarray(self.MEMS_IC, jnp.float32),
h=float(self.h),
Vin=jnp.asarray(Vin, jnp.float32),
theta=float(self.theta), # keyword-only & static
N=int(self.N), # keyword-only & static
sd=self.sd, # keyword-only & static
amp=float(self.amplification),
VDC=float(self.VDC),
stride=self.stride,
tau=self.tau_steps,
pos_buf=jnp.asarray(pos_buf, jnp.float32),
vel_buf=jnp.asarray(vel_buf, jnp.float32),
fb_gain=float(self.fb_gain),
buf_idx=jnp.asarray(buf_idx),
buf_cnt=jnp.asarray(buf_cnt)
)
if update_state:
self._set_state(mode, env_idx, pos_buf, vel_buf, buf_idx, buf_cnt)
if direct_fb:
# add a bias to neuron readings of u and an additional 1
# neuron vals + [u] + [1]
theta_vals = np.array(theta_vals).reshape(1,-1)
# neuron_vals = np.append(theta_vals,np.append(u,[[1]],axis=1), axis=1)
neuron_vals = np.concatenate([theta_vals, u, np.ones((1, 1))], axis=1)
return neuron_vals
if full_data:
return positions, theta_vals, Vin
else:
return np.array(theta_vals)
def _zero_reservoir(self, mode: str, env_idx=None):
state = self.res_state[mode]
if env_idx is None:
state["pos_bufs"] = state["pos_bufs"].at[:].set(0)
state["vel_bufs"] = state["vel_bufs"].at[:].set(0)
state["buf_idxs"] = state["buf_idxs"].at[:].set(0)
state["buf_cnts"] = state["buf_cnts"].at[:].set(0)
else:
env_idx = jnp.asarray(env_idx)
if env_idx.dtype == jnp.bool_:
env_idx = jnp.where(env_idx)[0]
state["pos_bufs"] = state["pos_bufs"].at[env_idx].set(0)
state["vel_bufs"] = state["vel_bufs"].at[env_idx].set(0)
state["buf_idxs"] = state["buf_idxs"].at[env_idx].set(0)
state["buf_cnts"] = state["buf_cnts"].at[env_idx].set(0)
self.res_state[mode] = state
def _convert_old_reservoir_state(self):
if hasattr(self, 'pos_buf'):
self.res_state = {
"train": {
"pos_bufs": jnp.expand_dims(jnp.asarray(self.pos_buf), axis=0),
"vel_bufs": jnp.expand_dims(jnp.asarray(self.vel_buf), axis=0),
"buf_idxs": jnp.array([self.buf_idx], dtype=jnp.int32),
"buf_cnts": jnp.array([self.buf_cnt], dtype=jnp.int32),
},
"val": {
"pos_bufs": jnp.expand_dims(jnp.asarray(self.pos_buf), axis=0),
"vel_bufs": jnp.expand_dims(jnp.asarray(self.vel_buf), axis=0),
"buf_idxs": jnp.array([self.buf_idx], dtype=jnp.int32),
"buf_cnts": jnp.array([self.buf_cnt], dtype=jnp.int32),
},
"test": {
"pos_bufs": jnp.expand_dims(jnp.asarray(self.pos_buf), axis=0),
"vel_bufs": jnp.expand_dims(jnp.asarray(self.vel_buf), axis=0),
"buf_idxs": jnp.array([self.buf_idx], dtype=jnp.int32),
"buf_cnts": jnp.array([self.buf_cnt], dtype=jnp.int32),
},
}
def _zero_reservoir_old(self):
self.pos_buf = jnp.full((self.buf_len,), 0.0, dtype=jnp.float32)
self.vel_buf = jnp.full((self.buf_len,), 0.0, dtype=jnp.float32)
self.buf_idx = jnp.int32(0)
self.buf_cnt = jnp.int32(0)
def _get_neuron_sat(self, neurons):
return np.sum(np.abs(neurons)>=0.89)*100/neurons.size
[docs]
@staticmethod
@partial(jax.jit, static_argnames=('n_steps','h','theta','N', 'sd', 'stride', 'tau', 'fb_gain', 'amp', 'VDC'))
def rk4(*, n_steps, pos_buf, vel_buf, h, Vin, theta, N, sd, amp, stride, tau, fb_gain, buf_idx, buf_cnt, VDC):
"""
Simulate and sample a DDE for a given set of step input voltages with RK4
and cubic Hermite interpolation for the delay term. The state of the DDE
is stored in ring buffers to allow for efficient and JAX correct sampling
and delayed feedback.
Args:
n_steps (int): total number of simulation steps.
pos_buf (jnp.ndarray): ring buffer to store past positions, length buf_len.
vel_buf (jnp.ndarray): ring buffer to store past velocities, length buf_len.
h (float): step size for the simulation.
Vin (jnp.ndarray): input voltage array, length n_steps.
theta (float): sample time for one reservoir reading [non-dimensional].
N (int): number of reservoir states.
sd (int): sample delay (in timesteps).
stride (int): sample every 'stride' with delay 'start'.
tau (float): delayed feedback gain parameter [non-dimensional].
fb_gain (float): feedback gain parameter [non-dimensional].
amp (float): amplitude of the input voltage [non-dimensional].
buf_idx (int): next write index (most recent sample is at buf_idx-1).
buf_cnt (int): number of valid samples in buffers (<= len(pos_buf)).
VDC (float): DC offset voltage [V].
Returns:
positions (jnp.ndarray): all positions from simulation.
sampled (jnp.ndarray): sampled positions (neuron values), length n_samples.
pbuf_out (jnp.ndarray): updated ring buffer of past positions.
vbuf_out (jnp.ndarray): updated ring buffer of past velocities.
buf_idx_out (int): updated next write index.
buf_count_out (int): updated number of valid samples in buffers.
"""
# time grid (for indices only)
tgrid = jnp.arange(n_steps, dtype=jnp.int32) # 0 ... n_steps-1
th = tgrid * h
# precompute Vin indices for k1 (t), k2/k3 (t+h/2), k4 (t+h)
idx1 = jnp.floor_divide(th, theta).astype(jnp.int32)
idx2 = jnp.floor_divide(th + 0.5*h, theta).astype(jnp.int32)
idx4 = jnp.floor_divide(th + h, theta).astype(jnp.int32)
# clamp indices to valid range
vmax = jnp.int32(Vin.shape[0] - 1)
idx1 = jnp.clip(idx1, 0, vmax)
idx2 = jnp.clip(idx2, 0, vmax)
idx4 = jnp.clip(idx4, 0, vmax)
# gather per-step voltages (all on device)
V1 = Vin[idx1] # [n_steps]
V2 = Vin[idx2]
V4 = Vin[idx4]
# init state
buf_len = pos_buf.shape[0]
last_idx = jnp.mod(buf_idx-1, buf_len)
y0f = jnp.asarray((pos_buf[last_idx],vel_buf[last_idx]), jnp.float32)
h32 = jnp.float32(h)
fb_gain = jnp.float32(fb_gain)
use_fb = tau > 0
tau = jnp.float32(tau)
# in discrete domain where h has been transformed into ints
# tau is h steps back of delay
# ring buffers
m = jnp.floor(tau).astype(jnp.int32)
init_pos_buf = jnp.asarray(pos_buf, jnp.float32) # jnp.full((buf_len,), y0f[0], dtype=jnp.float32) # prefill with initial pos
init_vel_buf = jnp.asarray(vel_buf, jnp.float32)
# carry.shape = (y, pos_ringbuf, i)
init_carry = (y0f, init_pos_buf, init_vel_buf, jnp.int32(0), jnp.int32(buf_idx), jnp.int32(buf_cnt))
def hermite_unit(x0, x1, v0, v1, unit):
# Cubic Hermite on unit interval, unit in [0,1]
# unit = (time - tau) - floor(time - tau) roughly
# H1 = 2*unit^3 - 3*unit^2 + 1
# H2 = -2*unit^3 + 3*unit^2
# H3 = unit^3 - 2*unit^2 + unit
# H4 = unit^3 - unit^2
# return H1*x0 + H2*x1 + H3*v0 + H4*v1
t2 = unit * unit
t3 = t2 * unit
H1 = 2.0*t3 - 3.0*t2 + 1.0
H2 = -2.0*t3 + 3.0*t2
H3 = t3 - 2.0*t2 + unit
H4 = t3 - t2
return H1*x0 + H2*x1 + H3*v0 + H4*v1
def sample_delay_hermite_ring(i, c, tau, pos_buf, vel_buf, buf_idx, buf_cnt, buf_len, m):
"""
Evaluate delayed position x(t_i + c - tau_steps) from a ring buffer using cubic Hermite.
pos_buf/vel_buf: ring buffers storing past samples at unit step spacing.
buf_idx: next write index (most recent sample is at buf_idx-1).
buf_cnt: number of valid samples in buffers (<= len(pos_buf)).
"""
# s = i + c - tau which is the s = timestep + rk4 offset - tau
s = (i.astype(jnp.float32) + jnp.float32(c)) - jnp.float32(tau)
j = jnp.floor(s).astype(jnp.int32) # left node
unit = s - jnp.floor(s) # fractional in [0,1)
# ex. i=0, c=0, tau=10, s=-10, j=-10, unit = -10-(-10) = 0
# ex. i=0, c=0.5, tau=10, s=-9.5, j=-9, unit = -9.5-(-10) = 0.5
# ex. i=0, c=1, tau=10, s=-9, j=-9, unit = -9-(-9) = 0
# ex. i=3, c=0.5, tau=7, s=-3.5 j=-4, unit = -3.5-(-4) = 0.5
# ex. i=7, c=0.5, tau=7, s=0.5, j=0, unit = 0.5-0 = 0.5
# Map logical index j to ring-buffer index. Latest logical index is i-1 at (buf_idx-1).
# So ring index for j is: buf_idx + j - i (mod buf_len). For j+1: +1.
j0 = (buf_idx + j - i) % buf_len
j1 = (buf_idx + j + 1 - i) % buf_len
x0 = pos_buf[j0]; x1 = pos_buf[j1]
v0 = vel_buf[j0]; v1 = vel_buf[j1] # velocity = dx/dt
# We need at least m+1 samples available to safely interpolate,
# where m = floor(tau_steps). This covers the worst case at c=0, i=0.
have_hist = buf_cnt >= (m + 1)
# Fallback (early warm-up): simple midpoint linear estimate
x_lin = 0.5 * (x0 + x1)
# x_h = hermite_unit(x0, x1, v0, v1, unit)
x_h = hermite_unit(x0, x1, v0*h32, v1*h32, unit)
return jax.lax.select(have_hist, x_h, x_lin)
def step(carry, Vs):
y, pbuf, vbuf, i, widx, cnt = carry
V1_t, V2_t, V4_t = Vs
if use_fb:
x_tau_c0 = sample_delay_hermite_ring(i, 0.0, tau, pbuf, vbuf, widx, cnt, buf_len, m)
x_tau_c05 = sample_delay_hermite_ring(i, 0.5, tau, pbuf, vbuf, widx, cnt, buf_len, m)
x_tau_c1 = sample_delay_hermite_ring(i, 1.0, tau, pbuf, vbuf, widx, cnt, buf_len, m)
fb_pos1 = fb_gain * x_tau_c0
fb_pos2 = fb_gain * x_tau_c05
fb_pos4 = fb_gain * x_tau_c1
else:
# feedback = 0 if tau=0
fb_pos1 = fb_pos2 = fb_pos4 = jnp.float32(0.0)
# apply feedback
V1_eff = V1_t + fb_pos1
V2_eff = V2_t + fb_pos2
V4_eff = V4_t + fb_pos4
k1 = h32 * Reservoir.diff_eq(y, V1_eff, amp, VDC)
k2 = h32 * Reservoir.diff_eq(y + 0.5 * k1, V2_eff, amp, VDC)
k3 = h32 * Reservoir.diff_eq(y + 0.5 * k2, V2_eff, amp, VDC)
k4 = h32 * Reservoir.diff_eq(y + k3, V4_eff, amp, VDC)
y_next = y + (k1 + 2*k2 + 2*k3 + k4) / 6.0
# stopper at pos >= 0.9
pos, vel = y_next[0], y_next[1]
clamp = pos >= 0.9
pos = jnp.where(clamp, jnp.float32(0.9), pos)
# vel = jnp.where(clamp, jnp.float32(0.0), vel)
vel = jnp.where(clamp & (vel > 0), jnp.float32(0.0), vel)
y_next = jnp.array([pos, vel], dtype=jnp.float32)
pbuf_next = pbuf.at[widx].set(y_next[0])
vbuf_next = vbuf.at[widx].set(y_next[1])
widx_next = jnp.mod(widx+1, buf_len)
cnt_next = jnp.minimum(cnt+1, jnp.int32(buf_len))
return (y_next, pbuf_next, vbuf_next, i+1, widx_next, cnt_next), y_next
# run scan
Vs_all = (V1, V2, V4) # tuple of [n_steps] arrays
(y_last, pbuf_out, vbuf_out, _, buf_idx_out, buf_count_out), Ys = jax.lax.scan(step, init_carry, Vs_all) # Ys: [n_steps, 2]
# f, init carry, xs
# .scan scans (iterates?) a function over an array while carrying along state
# sample every 'stride' with delay 'start'
positions = Ys[:, 0] # [n_steps]
sampled = positions[sd:sd+N*stride:stride]
# [n_samples]
buf_idx = jnp.mod(buf_idx, buf_len)
return (
positions,
sampled,
use_fb,
pbuf_out,
vbuf_out,
buf_idx_out,
buf_count_out
)
[docs]
@staticmethod
def diff_eq(y, V, amp, VDC):
"""
MEMS differential equation.
Args:
y (jnp.ndarray): Current state of the MEMS device [position, velocity].
V (jnp.ndarray): Input voltage at current time [Volts].
amp (jnp.float): Amplification parameter for the reservoir.
"""
zeta = jnp.float32(0.2)
eps = jnp.float32(8.85e-12)
k = jnp.float32(215.0)
A = jnp.float32(39.6e-6)
d = jnp.float32(42e-6)
VDC = jnp.float32(VDC)
Amp = jnp.float32(amp)
pos, vel = y[0], y[1]
den = 2.0 * k * (d**3) * jnp.maximum(1e-6, (1.0 - pos)**2)
acc = (eps * A * (Amp * V + VDC)**2) / den - 2.0 * zeta * vel - pos
return jnp.array([vel, acc], dtype=jnp.float32)
[docs]
class AnalyzeReservoir:
def __init__(self, reservoir):
self.reservoir = reservoir
[docs]
def plot_mask(self, *, together=True, separate=True, combined=True,
cmap_color='tab10', title_addition="", save=False, folder_path=None,
base_file_name="", figsize=(3.5,3), ax=None, show=True, color='royalblue',
partial_N=None):
if ax is None:
fig, ax = plt.figure(figsize=figsize, dpi=150)
else:
fig = ax.figure
if folder_path is None and save:
print('No folder provided, plot will not be saved.')
save = False
if folder_path is not None:
if not os.path.exists(folder_path):
os.makedirs(folder_path)
cmap = plt.get_cmap(cmap_color)
colors = [cmap(i) for i in range(self.reservoir.state_shape)]
if together:
for i in range(self.reservoir.state_shape):
y_axis = np.repeat(self.reservoir.mask[i, :], 50)
x_axis = np.linspace(0, self.reservoir.N-1, len(y_axis))
plt.plot(x_axis, y_axis, label=f'Mask {i+1}', color=colors[i])
plt.xlabel('Neuron')
plt.ylabel('Weight')
plt.title(f'Mask for Each State {title_addition}') if title_addition != "" else None
plt.legend()
plt.grid()
if save:
plt.savefig(f'{folder_path}/{base_file_name}_together_masks.png')
if separate:
for i in range(self.reservoir.state_shape):
plt.figure(figsize=figsize, dpi=150)
y_axis = np.repeat(self.reservoir.mask[i, :], 50)
x_axis = np.linspace(0, self.reservoir.N-1, len(y_axis))
plt.plot(x_axis, y_axis, label=f'Mask {i+1}', color=colors[i])
plt.plot(self.reservoir.mask[i, :], color='black')
plt.xlabel('Neuron')
plt.ylabel('Weight')
plt.title(f'Mask {i+1} {title_addition}') if title_addition != "" else None
plt.legend()
plt.grid()
if save:
plt.savefig(f'{folder_path}/{base_file_name}_mask_{i+1}.png')
if combined:
if partial_N is not None:
mask_partial = self.reservoir.mask[:, :partial_N]
else:
partial_N = self.reservoir.N
mask_partial = self.reservoir.mask
y_axis = np.repeat(mask_partial.sum(axis=0), 50)
x_axis = np.linspace(0, partial_N, len(y_axis))
ax.plot(x_axis, y_axis, label='Combined Mask', color=color, linewidth=1)
# ax.set_xlabel('Neuron')
ax.set_ylabel('Mask Weight', labelpad=0)
ax.set_xticks([0, 10, 20])
ax.set_xticklabels(['0', '10', '20',], fontsize=8)
ax.set_ylim(-1.1, 1.1)
ax.set_yticks([-1, -0.5, 0, 0.5, 1])
ax.set_xlim(0, partial_N)
ax.set_title('a) Combined Mask', pad=6, loc='left', fontsize=9)
# plt.legend()
ax.grid(True, linestyle='--', alpha=0.5)
# fig.tight_layout()
if save:
fig.savefig(f'{folder_path}/{base_file_name}_combined_mask.png')
if show:
plt.show()
[docs]
def plot_Vin(self, *, title_addition="", save=False, folder_path=None,
base_file_name="", figsize=(3.5,3), ax=None, show=True, color='black',
state=[0,0,0,0], Vin=None, partial_N=None, state2=None, color2='orange'):
if ax is None:
fig, ax = plt.figure(figsize=figsize, dpi=150)
else:
fig = ax.figure
if Vin is None:
Vin = self.reservoir.sim(obs=state, get_Vin=True)
if partial_N is not None:
Vin_partial = Vin[:partial_N]
else:
partial_N = self.reservoir.N
Vin_partial = Vin
y_axis = np.repeat(Vin_partial, 1000 // len(Vin_partial))
x_axis = np.linspace(0, partial_N, len(y_axis))
ax.plot(x_axis, y_axis, color=color, linewidth=1, linestyle='-')
ax.text(7, 0.4*360+20, r'$\mathbf{V}^{(0)}_{\text{in}}$', fontsize=8, ha='center', color=color)
if state2 is not None:
Vin2 = self.reservoir.sim(obs=state2, get_Vin=True)
if partial_N is not None:
Vin_partial = Vin2[:partial_N]
else:
partial_N = self.reservoir.N
Vin_partial = Vin
y_axis = np.repeat(Vin_partial, 1000 // len(Vin_partial))
x_axis = np.linspace(0, partial_N, len(y_axis))
ax.plot(x_axis, y_axis, color=color2, linewidth=1, linestyle='-', alpha=1)
ax.text(16, -0.3*360+20, r'$\mathbf{V}^{(1)}_{\text{in}}$', fontsize=8, ha='center', color=color2)
# ax.set_xlabel('Neuron')
ax.set_xticks([0, 10, 20])
ax.set_xticklabels(['0', '10', '20',], fontsize=8)
# ax.set_ylim(-1.1, 1.1)
ax.set_xlim(0, partial_N)
# ax.yaxis.tick_right()
# ax.yaxis.set_label_position('right')
# ax.set_yticklabels([])
ax.set_ylabel(r'Voltage $\left[ \text{V} \right]$ ', labelpad=0)
ax.set_title(r'b) Voltage Input', pad=6, loc='left', fontsize=9) # $\mathbf{V}_{\text{in}}$
# plt.legend()
ax.grid(True, linestyle='--', alpha=0.5)
# fig.tight_layout()
if save:
fig.savefig(f'{folder_path}/{base_file_name}_combined_mask.png')
if show:
plt.show()
[docs]
def plot_response(self, data, *, points=None, title=None, save=False,
folder_path=None, base_file_name="",
overlay_mask=False, Vin=None, show=True, figsize=(10, 5),
show_neuron_sat=False, neuron_sat=None, reps,
return_image=False, ax=None, color='royalblue', color2=None,
x_on=True, one_y=False, conf=False, Vin_color='green',
show_points=True, adjust_dict=None, xticks=None, xtick_labels=None):
if ax is None:
fig, ax = plt.subplots(figsize=figsize, dpi=150)
else:
fig = ax.figure
if folder_path is None and save:
print('No folder provided, plot will not be saved.')
save = False
if folder_path is not None:
if not os.path.exists(folder_path):
os.makedirs(folder_path)
x_axis = np.arange(len(data))/self.reservoir.stride
# x_axis = np.linspace(0, self.reservoir.N, len(data))
# plt.figure(figsize=figsize)
if color2 is not None:
# print(data.shape)
# print(data[:len(x_axis)//2].shape)
# print(data[len(x_axis)//2:].shape)
ax.plot(x_axis[:len(x_axis)//2], data[:len(x_axis)//2], color=color, linewidth=1, zorder=2)
ax.plot(x_axis[len(x_axis)//2:], data[len(x_axis)//2:], color=color2, linewidth=1, zorder=2)
else:
ax.plot(x_axis, data, color=color, label='Response', linewidth=1, zorder=2)
if points is not None and show_points:
# print(self.reservoir.sd/self.reservoir.stride)
x_axis = np.arange(len(points))+self.reservoir.sd/self.reservoir.stride
# print(points)
# x_axis = np.repeat()
# x_axis = np.linspace(0+self.reservoir.sd/self.reservoir.stride, self.reservoir.N, len(points))
ax.scatter(x_axis,points, color='black', marker='o', s=4, zorder=3, label='Neurons')
# print(len(points))
if overlay_mask:
y_axis = np.repeat(self.reservoir.mask.sum(axis=0), len(data) // len(self.reservoir.mask.sum(axis=0)))
x_axis = np.linspace(0, self.reservoir.N*reps, len(y_axis))
ax.plot(x_axis, y_axis, label='Combined Mask', color='red', linestyle='--', linewidth=1, zorder=2)
if Vin is not None:
# print(Vin)
y_axis = np.repeat(Vin, len(data) // len(Vin))
x_axis = np.linspace(0, self.reservoir.N*reps, len(y_axis))
ax.plot(x_axis, y_axis, label='Normed Voltage Input', color=Vin_color, linewidth=1, linestyle='-', zorder=0)
if show_neuron_sat:
ax.text(-1/32*self.reservoir.N, max(data)+.11, f'Neuron Saturation: {neuron_sat:.1f} %', fontsize=10, zorder=4)
if xticks is not None:
ax.set_xticks(xticks)
tick_labels = [str(int(round(tick*self.reservoir.theta, 0))) for tick in xticks]
ax.set_xticklabels(tick_labels, fontsize=8)
if x_on:
ax.set_xlabel(r'Non-dimensional Time, $\tilde{t}$', fontsize=8, labelpad=2)
# else:
# ax.set_xticklabels([], fontsize=8)
if not one_y:
ax.set_ylabel(r'MEMS Response, $\tilde{x}$', fontsize=8, labelpad=0)
ax.set_ylim([-0.5, 1.1])
ax.set_yticks([-0.5, 0, 0.5, 1])
ax.set_yticklabels(['-0.5', '0', '0.5', '1'], fontsize=8)
if title is not None:
ax.set_title(title, pad=6, loc='left', fontsize=9)
# plt.legend()
# plt.tight_layout(pad=.5)
# ax.axhline(y=0.9, color='red', linestyle='--', zorder=1, linewidth=1)
if conf:
ax.axhline(y=0.9, xmin=0, xmax=.05, color='red', linestyle='--', zorder=1, linewidth=1)
ax.axhline(y=0.9, xmin=.11, xmax=1, color='red', linestyle='--', zorder=1, linewidth=1)
ax.text(1.35, 0.89, r'$x_s$', fontsize=8, zorder=1, color='red',
ha='center', va='center', fontweight='bold')
vline_color = "#888888"
x_color = "#444444"
ax.axvline(x=20, color=vline_color, linestyle='--', zorder=1, linewidth=1)
ax.text(1.3, -0.4, r'$\mathbf{x}_0$', fontsize=8, zorder=1, color=x_color,ha='center')
ax.text(21.3, -0.4, r'$\mathbf{x}_1$', fontsize=8, zorder=1, color=x_color,ha='center')
# ax.annotate('90% Neuron Saturation', xy=(self.reservoir.N/2, 0.9), xytext=(self.reservoir.N/2+2, 0.9), ha='center', va='center', fontsize=7, color='red')
# plt.text(0, 0.77, r'$x_s$', fontsize=10, zorder=1, color='red')
ax.grid(True, linestyle='--', alpha=0.5)
self._subplots_adjust(fig, adjust_dict=adjust_dict)
if save:
plt.savefig(f'{folder_path}/{base_file_name}_sim_response.png')
if show:
plt.show()
if return_image:
buf = BytesIO()
plt.savefig(buf, format='png', dpi=75)
buf.seek(0)
image = Image.open(buf)
plt.close()
return image
else:
return None
# plt.figure(figsize=figsize)
# plt.plot(data)
# plt.show()
[docs]
def sim_response(self, *, state=None, seed=None, overlay_mask=False,
overlay_Vin=False, overlay_actual_Vin=True,
title=None, save=False, folder_path=None, base_file_name="",
show=True, figsize=(10, 5), show_neuron_sat=False, reps=1,
return_image=False, states=None, ax=None, color='royalblue', color2='royalblue',
x_on=True, one_y=False, Vin_color='green', show_points=True,
adjust_dict=None, xticks=None, xtick_labels=None):
# Vin = np.array([0])
if state is None:
rng = np.random.RandomState(seed)
state = rng.uniform(0, 1, self.reservoir.state_shape)
# print(f'Random state: {state}')
all_positions = np.array([])
all_theta_vals = np.array([])
all_Vin = np.array([])
for i in range(reps):
if states is not None:
state = states[i]
elif states is None:
if i>0:
state = rng.uniform(0, 1, self.reservoir.state_shape)
positions, theta_vals, Vin = self.reservoir.sim(obs=state, full_data=True)
all_positions = np.append(all_positions, positions)
all_theta_vals = np.append(all_theta_vals, theta_vals)
all_Vin = np.append(all_Vin, Vin)
# print(all_theta_vals)
positions = np.reshape(all_positions, (-1, ))
theta_vals = np.reshape(all_theta_vals, (-1,))
Vin = np.reshape(all_Vin, (-1,))
# print(len(theta_vals))
# positions, theta_vals, Vin = self.reservoir.sim(obs=state, full_data=True)
neuron_sat = self.reservoir._get_neuron_sat(theta_vals)
# print(np.mean(positions))
# print(np.std(positions))
if not overlay_Vin:
Vin = None
if overlay_actual_Vin and overlay_Vin:
Vin = abs(Vin*self.reservoir.amplification + self.reservoir.VDC)
Vin = Vin / max(Vin)
image = self.plot_response(positions, points=theta_vals, overlay_mask=overlay_mask,
Vin=Vin, title=title, save=save, folder_path=folder_path,
base_file_name=base_file_name, show=show, figsize=figsize,
show_neuron_sat=show_neuron_sat, neuron_sat=neuron_sat,
reps=reps, return_image=return_image, ax=ax, color=color, color2=color2,
x_on=x_on, one_y=one_y, Vin_color=Vin_color, show_points=show_points,
adjust_dict=adjust_dict, xticks=xticks, xtick_labels=xtick_labels)
return image
def _subplots_adjust(self, fig, adjust_dict=None, *, left=0.13, right=0.976, top=0.936, bottom=0.107, wspace=0.438, hspace=0.35):
for key, value in adjust_dict.items():
if value is not None:
if key == 'left':
left = value
elif key == 'right':
right = value
elif key == 'top':
top = value
elif key == 'bottom':
bottom = value
elif key == 'wspace':
wspace = value
elif key == 'hspace':
hspace = value
fig.subplots_adjust(left=left, right=right, top=top, bottom=bottom, wspace=wspace, hspace=hspace)
[docs]
def reservoir_subplots(self, *, figsize=(3.5, 2.5), dpi=150,states=None,
show=True, save=True, folder_path=None, file_name=None):
plt.rcParams['font.sans-serif'] = "Arial"
plt.rcParams['font.family'] = "sans-serif"
plt.rcParams.update({
"font.size": 8,
"axes.labelsize": 9,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"legend.fontsize": 8,
})
fig = plt.figure(figsize=figsize, dpi=dpi)
gs = fig.add_gridspec(2, 2, width_ratios=[1, 1], height_ratios=[1,1])
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[1, 0:])
self.plot_mask(ax=ax1, save=False, show=False, together=False, separate=False,
combined=True, color='black', partial_N=20)
color='#0072B2'
color2='#D55E00'
self.plot_Vin(ax=ax2, save=False, show=False, color=color, color2=color2,state=[0,0,0,0],
partial_N=20, state2=[-.1, -0.5, -0.15, -0.5])
self.sim_response(ax=ax3, save=False, show=False, states=[[0,0,0,0], [-0.1, -0.5, -0.15, -0.5]],
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2)
fig.subplots_adjust(left=0.13,
right=0.976,
top=0.936,
bottom=0.107,
wspace=0.438,
hspace=0.35,
)
if save:
fig.savefig(f'{folder_path}/{file_name}.png', dpi=300)
if show:
plt.show()
return fig
[docs]
def reservoir_response_subplots_3(self, *, figsize=(3.5, 2.5), dpi=150,states=None,
show=True, save=True, folder_path=None, file_name=None):
plt.rcParams['font.sans-serif'] = "Arial"
plt.rcParams['font.family'] = "sans-serif"
plt.rcParams.update({
"font.size": 8,
"axes.labelsize": 9,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"legend.fontsize": 8,
})
fig = plt.figure(figsize=figsize, dpi=dpi)
gs = fig.add_gridspec(3, 1, height_ratios=[1,1,1])
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[2, 0])
color='#0072B2'
color2='#D55E00'
titles = [r'a) $\theta$ = 0.1', r'b) $\theta$ = 1', r'c) $\theta$ = 2$\pi$']
self.reservoir = Reservoir(N=20,
theta=0.1, #*np.pi,
state_shape=4,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [4.8*2,3+4,2*0.418,2+4], #changed
norm_offset= [4.8,4,0.5,4],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax1, save=False, show=False, states=[[0,0,0,0], [-0.1, -0.5, -0.15, -0.5]],
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[0], one_y=True)
self.reservoir = Reservoir(N=20,
theta=1, #*np.pi,
state_shape=4,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [4.8*2,3+4,2*0.418,2+4], #changed
norm_offset= [4.8,4,0.5,4],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax2, save=False, show=False, states=[[0,0,0,0], [-0.1, -0.5, -0.15, -0.5]],
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[1], one_y=True)
self.reservoir = Reservoir(N=20,
theta=2*np.pi,
state_shape=4,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [4.8*2,3+4,2*0.418,2+4], #changed
norm_offset= [4.8,4,0.5,4],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax3, save=False, show=False, states=[[0,0,0,0], [-0.1, -0.5, -0.15, -0.5]],
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2,
title=titles[2], one_y=True)
fig.subplots_adjust(left=0.13,
right=0.976,
top=0.936,
bottom=0.107,
wspace=0.438,
hspace=0.552,
)
fig.text(0.01, 0.5, r'MEMS Response, $\tilde{x}$',
va='center', rotation='vertical')
if save:
fig.savefig(f'{folder_path}/{file_name}.png', dpi=300)
if show:
plt.show()
return fig
[docs]
def reservoir_response_subplots_4(self, *, figsize=(3.5, 2.5), dpi=150,states=None,
show=True, save=True, folder_path=None, file_name=None):
plt.rcParams['font.sans-serif'] = "Arial"
plt.rcParams['font.family'] = "sans-serif"
plt.rcParams.update({
"font.size": 8,
"axes.labelsize": 9,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"legend.fontsize": 8,
})
fig = plt.figure(figsize=figsize, dpi=dpi)
gs = fig.add_gridspec(4, 1, height_ratios=[1,1,1,1])
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[2, 0])
ax4 = fig.add_subplot(gs[3, 0])
color='#0072B2'
color2='#D55E00'
titles = [r'a) $\theta$ = 0.1', r'b) $\theta$ = 0.5', r'c) $\theta$ = 1', r'd) $\theta$ = 2$\pi$']
self.reservoir = Reservoir(N=20,
theta=0.1, #*np.pi,
state_shape=4,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [4.8*2,3+4,2*0.418,2+4], #changed
norm_offset= [4.8,4,0.5,4],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax1, save=False, show=False, states=[[0,0,0,0], [-0.1, -0.5, -0.15, -0.5]],
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[0], one_y=True)
self.reservoir = Reservoir(N=20,
theta=0.5, #*np.pi,
state_shape=4,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [4.8*2,3+4,2*0.418,2+4], #changed
norm_offset= [4.8,4,0.5,4],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax2, save=False, show=False, states=[[0,0,0,0], [-0.1, -0.5, -0.15, -0.5]],
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[1], one_y=True)
self.reservoir = Reservoir(N=20,
theta=1, #*np.pi,
state_shape=4,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [4.8*2,3+4,2*0.418,2+4], #changed
norm_offset= [4.8,4,0.5,4],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax3, save=False, show=False, states=[[0,0,0,0], [-0.1, -0.5, -0.15, -0.5]],
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2,
title=titles[2], x_on=False, one_y=True)
self.reservoir = Reservoir(N=20,
theta=2*np.pi,
state_shape=4,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [4.8*2,3+4,2*0.418,2+4], #changed
norm_offset= [4.8,4,0.5,4],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax4, save=False, show=False, states=[[0,0,0,0], [-0.1, -0.5, -0.15, -0.5]],
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2,
title=titles[3], one_y=True)
fig.subplots_adjust(left=0.13,
right=0.976,
top=0.954,
bottom=0.087,
wspace=0.438,
hspace=0.624,
)
fig.text(0.01, 0.5, r'MEMS Response, $\tilde{x}$',
va='center', rotation='vertical')
if save:
fig.savefig(f'{folder_path}/{file_name}.png', dpi=300)
if show:
plt.show()
return fig
[docs]
def reservoir_response_subplots_2_mc(self, *, figsize=(3.5, 2.5), dpi=150,states=None,
show=True, save=True, folder_path=None, file_name=None):
plt.rcParams['font.sans-serif'] = "Arial"
plt.rcParams['font.family'] = "sans-serif"
plt.rcParams.update({
"font.size": 8,
"axes.labelsize": 9,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"legend.fontsize": 8,
})
fig = plt.figure(figsize=figsize, dpi=dpi)
gs = fig.add_gridspec(2, 1, height_ratios=[1,1])
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[1, 0])
color='#0072B2'
color2='#D55E00'
titles = [r'a) $\tau$ = 0', r'b) $\tau$ = 20']
states = [[0,0], [-0.2, -0.05]]
self.reservoir = Reservoir(N=20,
theta=1, #*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax1, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[0])
self.reservoir = Reservoir(N=20,
theta=1, #*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=20,
fb_gain=0.5,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax2, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2,
title=titles[1])
fig.subplots_adjust(left=0.13,
right=0.976,
top=0.933,
bottom=0.12,
wspace=0.438,
hspace=0.3,
)
if save:
fig.savefig(f'{folder_path}/{file_name}.png', dpi=300)
if show:
plt.show()
return fig
[docs]
def reservoir_response_subplots_3_mc(self, *, figsize=(3.5, 2.5), dpi=150,states=None,
show=True, save=True, folder_path=None, file_name=None):
plt.rcParams['font.sans-serif'] = "Arial"
plt.rcParams['font.family'] = "sans-serif"
plt.rcParams.update({
"font.size": 8,
"axes.labelsize": 9,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"legend.fontsize": 8,
})
fig = plt.figure(figsize=figsize, dpi=dpi)
gs = fig.add_gridspec(3, 1, height_ratios=[1,1,1])
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[2, 0])
color='#0072B2'
color2='#D55E00'
titles = [r'a) $\theta$ = 0.1', r'b) $\theta$ = 0.5', r'c) $\theta$ = 1']
states = [[0,0], [-0.2, -0.05]]
self.reservoir = Reservoir(N=20,
theta=0.1, #*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax1, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[0],one_y=True)
self.reservoir = Reservoir(N=20,
theta=0.5, #*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax2, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[1],one_y=True)
self.reservoir = Reservoir(N=20,
theta=1, #*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=0,
fb_gain=0,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax3, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2,
title=titles[2],one_y=True)
fig.subplots_adjust(left=0.13,
right=0.976,
top=0.936,
bottom=0.107,
wspace=0.438,
hspace=0.35,
)
fig.text(0.01, 0.5, r'MEMS Response, $\tilde{x}$',
va='center', rotation='vertical')
if save:
fig.savefig(f'{folder_path}/{file_name}.png', dpi=300)
if show:
plt.show()
return fig
[docs]
def reservoir_response_subplots_4_mc(self, *, figsize=(3.5, 2.5), dpi=150,states=None,
show=True, save=True, folder_path=None, file_name=None):
plt.rcParams['font.sans-serif'] = "Arial"
plt.rcParams['font.family'] = "sans-serif"
plt.rcParams.update({
"font.size": 8,
"axes.labelsize": 9,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"legend.fontsize": 8,
})
fig = plt.figure(figsize=figsize, dpi=dpi)
gs = fig.add_gridspec(4, 1, height_ratios=[1,1,1,1])
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[2, 0])
ax4 = fig.add_subplot(gs[3, 0])
color='#0072B2'
color2='#D55E00'
titles = [r'a) $\theta$ = 0.1', r'b) $\theta$ = 0.5', r'c) $\theta$ = 1',r'd) $\theta$ = 2$\pi$']
states = [[0,0], [-0.2, -0.05]]
self.reservoir = Reservoir(N=20,
theta=0.1, #*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=20,
fb_gain=0.5,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax1, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[0],one_y=True)
self.reservoir = Reservoir(N=20,
theta=0.5, #*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=20,
fb_gain=0.5,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax2, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2, x_on=False,
title=titles[1],one_y=True)
self.reservoir = Reservoir(N=20,
theta=1, #*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=20,
fb_gain=0.5,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax3, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2,
title=titles[2], x_on=False,one_y=True)
self.reservoir = Reservoir(N=20,
theta=2*np.pi,
state_shape=2,
mask_seed=2,
h=0.02,
input_connectivity=0.2,
normalize_mask=True,
# norm_factor= [4.8*2,3+4,2*0.418,2+4],
# norm_offset=[4.8,4,0.5,4],
norm_factor= [2*1.2, 2*0.07],
norm_offset= [1.2, 0.07],
amp=360,
VDC=20,
tau=20,
fb_gain=0.5,
sd=0,
# load=True,
# mask=mask,
)
self.sim_response(ax=ax4, save=False, show=False, states=states,
overlay_mask=False, overlay_Vin=False, overlay_actual_Vin=False,
show_neuron_sat=False, reps=2,color=color,color2=color2,
title=titles[3],one_y=True)
fig.subplots_adjust(left=0.13,
right=0.976,
top=0.95,
bottom=0.087,
wspace=0.438,
hspace=0.68,
)
fig.text(0.01, 0.5, r'MEMS Response, $\tilde{x}$',
va='center', rotation='vertical')
if save:
fig.savefig(f'{folder_path}/{file_name}.png', dpi=300)
if show:
plt.show()
return fig
[docs]
def find_ss_amp(self, positions, delay):
from scipy.signal import find_peaks
# print(f'positions shape: {positions.shape}')
x = positions[delay:]
# print(f'x: {x}')
x_np = np.asarray(x)
peaks, _ = find_peaks(x_np)
troughs, _ = find_peaks(-x_np)
# print(f'peaks: {peaks}')
# print(f'troughs: {troughs}')
unique_peaks = np.unique(peaks)
unique_troughs = np.unique(troughs)
amp_peaks = x_np[unique_peaks]
amp_troughs = x_np[unique_troughs]
# print(amp_peaks)
# print(amp_troughs)
amp = np.mean(amp_peaks) - np.mean(amp_troughs)
# print(f'ss_amp: {amp}')
return amp
[docs]
def bode_plot(self, *, N=80, amplification=5, VAC=1, delay=35,
folder_path, file_name, show=True, save=False, fig_size=(5, 3.5)):
import matplotlib.ticker as ticker
import time
wns = jnp.logspace(-1, 1, 300)
freqs = wns/(2*np.pi)
ss_amps = []
for freq in freqs:
# ensure the resolution is sufficient
if freq >= 1/2/np.pi: # decrease h when increasing freq
steps = 5000
h = 1/(freq)/steps
print(f'freq={freq}, h={h}, steps={steps}')
else: # increase steps when decreasing freq
h = 0.002
steps = int(1/(freq)/h)
h = 1/(freq)/steps
print(f'freq={freq}, h={h}, steps={steps}')
# theta = np.ceil(steps*h)
theta = h
delay_steps = int(.75*steps)
Vin = self.create_sine_input(steps, VAC, N)
# print(f'Vin: {Vin.shape}')
reservoir = Reservoir(N=N, h=h, theta=theta, amp=amplification, n_steps_needed=N*steps)
positions, theta_vals = reservoir.sim(Vin=Vin, full_data=True)[:2]
# plt.figure()
# plt.plot(positions, color='black')
# # plt.plot(Vin, color='blue', alpha=0.5)
# plt.show()
# # time.sleep(1)
# plt.close()
# print(f'positions: {positions.shape}')
# print(f'theta_vals: {theta_vals.shape}')
ss_amp = self.find_ss_amp(positions, delay_steps)
print(f'freq: {freq:.4f}, ss_amp: {ss_amp:.4f}, delay_steps: {delay_steps}')
ss_amps.append([freq, ss_amp])
ss_amps = np.array(ss_amps)
plt.rcParams['font.sans-serif'] = "Arial"
plt.rcParams['font.family'] = "sans-serif"
plt.rcParams.update({
"font.size": 8,
"axes.labelsize": 9,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"legend.fontsize": 8,
})
plt.figure(figsize=fig_size, dpi=150)
plt.plot(2*np.pi*ss_amps[:,0], ss_amps[:,1], color='black')
plt.xlabel(r'Frequency Ratio, $\frac{\omega}{\omega_n}$')
plt.ylabel(r'Steady State Amplitude, $\tilde{x}$')
# plt.title(f'SS Amp (exact) Bode Plot, VAC={amplification*VAC}, VDC=20')
plt.xscale('log')
# plt.gca().xaxis.set_major_formatter(ticker.LogFormatterMathtext(base=10.0, labelOnlyBase=True))
plt.xticks([0.1, 1, 10])
plt.xlim([0.1, 10])
plt.grid(True, linestyle='--', alpha=0.5)
plt.text(.2, .009, f'Quasi-static', fontsize=8, color='black')
plt.text(2.5, .009, f'Inertia-dominated', fontsize=8, color='black')
plt.tight_layout()
plt.savefig(f'{folder_path}/{file_name}.png', dpi=300)
plt.show()