Skip to content

Commit

Permalink
Reparameterized, better design phase
Browse files Browse the repository at this point in the history 
I reparameterized the SYN stage in terms of spontaneous and saturation rates for the classes and few other things, and implemented the design process to account for everything to get those to come out right, all in terms of variables that show up on the diagram that's an augmentation of the v2 two_cap IHC diagram (for documentation, in progress).
Still need to move some of the numbers into SYN_params.
  • Loading branch information
@dicklyon
dicklyon committed Sep 19, 2024
1 parent 6975bde commit 17f5564
Show file tree
Hide file tree
Showing 5 changed files with 104 additions and 81 deletions.
81 changes: 54 additions & 27 deletions matlab/CARFAC_Design.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,8 +156,7 @@
'fs', fs, ...
'n_classes', n_classes, ...
'tau_lpf', 0.000080, ...
'out_rate', 0.02, ... % Depletion can be quick (few ms).
'recovery', 0.0001); % Recovery tau about 1000 to 10,000 sample times?
'reservoir_tau', 0.020);
end

% first figure out how many filter stages (PZFC/CARFAC channels):
Expand Down Expand Up @@ -579,43 +578,71 @@

n_ch = length(pole_freqs);
n_cl = SYN_params.n_classes;
col_n_ch_ones = ones(n_ch, 1);

v_width = 0.05;
v_widths = v_width * ones(1, n_cl);
max_rate = 2500; % % Instantaneous max at onset, per Kiang figure 5.18.
max_rates = max_rate * ones(1, n_cl);
sat_rates = 200;
switch n_cl % Just enough to test that both 2 and 3 can work.
case 2
offsets = [3, 6];
agc_weights = (fs/1000) * [0.125; 0.6];
res_inits = [0.2, 0.6];
rest_output = 0.8 * (fs/1000) * 0.016; % Subject off to get agc_in near 0 in quiet.
n_fibers = [50, 60]
v_half = v_widths .* [3, 6;]
spont_rates = [50, 5];
n_fibers = [50, 60];
agc_weights = [1.2, 2.4]'; % a column of weights
% Pick a sat level of reservoir state (move into params?):
w1 = 0.125 * [1, 1]; % How low the res gets in saturation.
case 3
offsets = [3, 5, 7];
agc_weights = (fs/1000) * [0.1; 0.5; 2];
res_inits = [0.13, 0.55, 0.9]
rest_output = (fs/1000) * 0.014;
spont_rates = [50, 6, 1];
n_fibers = [50, 35, 25];
v_half = v_widths .* [3, 5, 7]; % Compute from sponts instead?
agc_weights = [1.2, 1.2, 1.2]'; % a column of weights
% Pick a sat level of reservoir state (move into params?):
w1 = 0.2 * [1, 1, 1]; % How low the res gets in saturation.
otherwise
error('unimplemented n_classes in in SYN_params in CARFAC_DesignSynapses');
end
v_width = 0.02;
% Plural var names indicate values for the different rate classes.
v_widths = v_width * ones(1, n_cl);

% Copy stuff from params to coeffs and design a few things.
% Do some design. First, gains to get sat_rate when sigmoid is 1, which
% involves the reservoir steady-state solution.
% Assume the sigmoid is switching between 0 and 1 at 50% duty cycle, so
% normalized value is 0.5.
% Most of these are not per-channel, but could be expanded that way
% later if desired.

% Mean sat prob of spike per sample per neuron, likely same for all
% classes.
% Use names 1 for sat and 0 for spont in some of these.
p1 = sat_rates / fs;
p0 = spont_rates / fs;

q1 = 1 - w1;
s1 = 0.5; % average sigmoid output at saturation
z1 = s1*w1;
% solve q1 = a1*z1 for gain coeff a1:
a1 = q1 ./ z1;
% solve p1 = a2*z1 for gain coeff a2:
a2 = p1 ./ z1;

% Now work out how to get the desired spont.
z0 = p0 ./ a2;
q0 = z0 .* a1;
w0 = 1 - q0;
s0 = z0 ./ w0;
% Solve for (negative) sigmoid midpoint offset that gets s0 right.
offsets = log((1 - s0)./s0);

spont_p = a2 .* w0 .* s0; % should match p0; check it; yes it does.

% Copy stuff needed at run time into coeffs.
SYN_coeffs = struct( ...
'n_ch', n_ch, ...
'n_classes', n_cl, ...
'max_probs', max_rates / SYN_params.fs, ...
'n_fibers', col_n_ch_ones * n_fibers, ... % Synaptopathy comes in here; channelize it, too.
'v_width', v_width, ...
'v_half', col_n_ch_ones * v_half, ... % Same units as v_width and v_recep.
'out_rate', 0.1, ... % Depletion can be quick (few ms).
'recovery', 1e-3, ... % Recovery tau about 1000 sample times. Or 10X this?
'n_fibers', ones(n_ch,1) * n_fibers, ... % Synaptopathy comes in here.
'v_widths', v_widths, ...
'v_halfs', offsets .* v_widths, ... % Same units as v_recep and v_widths.
'a1', a1, ... % Feedback gain
'a2', a2, ... % Output gain
'agc_weights', agc_weights, ... % try to make syn_out resemble ihc_out to go to agc_in.
'rest_output', rest_output, ...
'res_inits', res_inits, ...
'spont_p', spont_p, ...
'res_lpf_inits', q0, ...
'res_coeff', 1 - exp(-1/(SYN_params.reservoir_tau * fs)), ...
'lpf_coeff', 1 - exp(-1/(SYN_params.tau_lpf * fs)));

7 changes: 4 additions & 3 deletions matlab/CARFAC_IHC_Step.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
% See the License for the specific language governing permissions and
% limitations under the License.

function [ihc_out, state, receptor_potential] = CARFAC_IHC_Step(bm_out, coeffs, state);
function [ihc_out, state, v_recep] = CARFAC_IHC_Step(bm_out, coeffs, state);
% function [ihc_out, state] = CARFAC_IHC_Step(bm_out, coeffs, state);
%
% One sample-time update of inner-hair-cell (IHC) model, including the
Expand All @@ -25,9 +25,9 @@
% receptor_potential output will be empty except in two_cap mode.
% Use it as input to CARFAC_SYN_Step to model synapses to get firing rates.

v_recep = []; % For cases other than two_cap and do_syn it's not used.
if coeffs.just_hwr
ihc_out = min(2, max(0, bm_out)); % limit it for stability
receptor_potential = [];
else
conductance = CARFAC_Detect(bm_out); % rectifying nonlinearity

Expand All @@ -44,7 +44,6 @@
state.lpf2_state = state.lpf2_state + coeffs.lpf_coeff * ...
(state.lpf1_state - state.lpf2_state);
ihc_out = state.lpf2_state - coeffs.rest_output;
receptor_potential = [];
else
% Change to 2-cap version mediated by receptor potential at cap1:
% Geisler book fig 8.4 suggests 40 to 800 Hz corner.
Expand All @@ -68,6 +67,8 @@
state.lpf1_state = state.lpf1_state + coeffs.lpf_coeff * ...
(ihc_out - state.lpf1_state);
ihc_out = state.lpf1_state - coeffs.rest_output;
% Return a modified receptor potential that's zero at rest, for SYN.
v_recep = coeffs.rest_cap1 - state.cap1_voltage;
end
end

Expand Down
5 changes: 2 additions & 3 deletions matlab/CARFAC_Init.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,6 @@
n_ch = coeffs.n_ch;
n_cl = coeffs.n_classes;
state = struct( ...
'reservoirs', ones(n_ch, 1) * coeffs.res_inits, ... % 1 means full.
'lpf_state', zeros(n_ch, n_cl) );
% Need to make different rest fullnesses for the different classes!!!
'reservoirs', ones(n_ch, 1) * coeffs.res_lpf_inits, ... % 0 full, 1 empty.
'lpf_state', ones(n_ch, 1) * coeffs.spont_p);

27 changes: 14 additions & 13 deletions matlab/CARFAC_SYN_Step.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,23 +3,24 @@
% returning instantaneous spike rates per class, for a group of neurons
% associated with the CF channel, including reductions due to synaptopathy.

% Normalized offset position into neurotransmitter release sigmoid.
x = (v_recep - coeffs.v_halfs) ./ coeffs.v_widths ;

% This sigmoidal release_rates is relative to a max of 1, usually way lower.
release_rates = state.reservoirs ./ ...
(1 + exp(-(v_recep - coeffs.v_half)./coeffs.v_width));
s = 1 ./ (1 + exp(-x));
release_rates = (1 - state.reservoirs) .* s; % z = w*s

% Smooth once with LPF (receptor potential was already smooth):
state.lpf_state = state.lpf_state + coeffs.lpf_coeff * ...
(release_rates - state.lpf_state);
release_rates = state.lpf_state;

% Deplete reservoirs some, independent of synaptopathy numbers (we assume).
state.reservoirs = state.reservoirs - coeffs.out_rate .* release_rates;
% Refill reservoirs a little toward 1.
state.reservoirs = state.reservoirs + coeffs.recovery .* (1 - state.reservoirs);

% Includes number of effective neurons per channel here, and interval T,
(coeffs.a2 .* release_rates - state.lpf_state); % this is firing probs.
firing_probs = state.lpf_state; % Poisson rate per neuron per sample.
% Include number of effective neurons per channel here, and interval T;
% so the rates (instantaneous action potentials per second) can be huge.
firings = coeffs.n_fibers .* coeffs.max_probs .* release_rates;
firings = coeffs.n_fibers .* firing_probs;

% Feedback, to update reservoir state q for next time.
state.reservoirs = state.reservoirs + coeffs.res_coeff .* ...
(coeffs.a1 .* release_rates - state.reservoirs);

% Make an output that resembles ihc_out, to go to agc_in.
syn_out = firings * coeffs.agc_weights - coeffs.rest_output;
syn_out = (coeffs.n_fibers .* (firing_probs - coeffs.spont_p)) * coeffs.agc_weights;
65 changes: 30 additions & 35 deletions matlab/CARFAC_Test.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,11 +157,11 @@
syn_state = CF.ears(1).SYN_state;
end
for k = 1:length(x_in)
[ihc_out, ihc_state, receptor_potential] = CARFAC_IHC_Step( ...
[ihc_out, ihc_state, v_recep] = CARFAC_IHC_Step( ...
x_in(k), CF.ears(1).IHC_coeffs, ihc_state);
if CF.do_syn % ignore ihc_out and use receptor_potential.
[syn_out, firings, syn_state] = CARFAC_SYN_Step( ...
receptor_potential, CF.ears(1).SYN_coeffs, syn_state);
v_recep, CF.ears(1).SYN_coeffs, syn_state);
% This can go a little negative; should be zero at rest.
neuro_output(k) = syn_out(1);
class_firings(k, :) = firings(1, :);
Expand Down Expand Up @@ -237,10 +237,8 @@
expected_maxes = expected_results(:, 1)';
expected_acs = expected_results(:, 2)';

fprintf(1, 'Golden data for Matlab test_IHC:\n');
fprintf(1, 'Golden data for Matlab and Python test_IHC:\n');
fprintf(1, ' [%f, %f];\n', [blip_maxes; blip_ac])
fprintf(1, 'Golden data for Python test_IHC:\n');
fprintf(1, ' [%f, %f],\n', [blip_maxes; blip_ac])

for i = 1:num_blips
if abs(expected_maxes(i) - blip_maxes(i)) > expected_maxes(i)/1e6
Expand Down Expand Up @@ -304,10 +302,8 @@
expected_maxes = expected_results(:, 1)';
expected_acs = expected_results(:, 2)';

fprintf(1, 'Golden data for Matlab test_IHC2:\n');
fprintf(1, 'Golden data for Matlab and Python test_IHC2:\n');
fprintf(1, ' [%f, %f];\n', [blip_maxes; blip_ac])
fprintf(1, 'Golden data for Python test_IHC2:\n');
fprintf(1, ' [%f, %f],\n', [blip_maxes; blip_ac])

for i = 1:num_blips
if abs(expected_maxes(i) - blip_maxes(i)) > expected_maxes(i)/1e6
Expand Down Expand Up @@ -340,21 +336,21 @@
switch test_freqs(k)
case 300
expected_results = [ ...
[0.747350, 177.802226];
[1.831880, 329.698994];
[3.635296, 517.895480];
[5.666424, 656.305962];
[7.522566, 721.465796];
[9.038596, 735.407690];
[0.582030, 101.529701];
[1.879711, 266.366280];
[3.399757, 461.743150];
[3.911907, 527.148854];
[3.915707, 511.042054];
[3.926793, 493.685648];
];
case 3000
expected_results = [ ...
[0.198211, 29.277988];
[0.523911, 55.506573];
[1.201355, 97.977400];
[2.264493, 146.819176];
[3.465151, 187.067551];
[4.500863, 212.210966];
[0.092435, 13.742136];
[0.342293, 40.817636];
[1.044103, 96.671170];
[1.952015, 148.463834];
[2.701081, 167.405953];
[3.071866, 153.483845];
];
otherwise
fprintf(1, 'No test_results for %f Hz in test_IHC.\n', ...
Expand All @@ -371,19 +367,19 @@
expected_maxes = expected_results(:, 1)';
expected_acs = expected_results(:, 2)';

fprintf(1, 'Golden data for Matlab test_IHC3:\n');
fprintf(1, 'Golden data for Matlab and Python test_IHC3:\n');
fprintf(1, ' [%f, %f];\n', [blip_maxes; blip_ac])
fprintf(1, 'Golden data for Python test_IHC3:\n');
fprintf(1, ' [%f, %f],\n', [blip_maxes; blip_ac])

for i = 1:num_blips
if abs(expected_maxes(i) - blip_maxes(i)) > expected_maxes(i)/1e6
% Lowered precision from 1e6 to 0.5e6 as the values are lower,
% and don't quite quite enough digits printed in the default format.
if abs(expected_maxes(i) - blip_maxes(i)) > expected_maxes(i)/0.5e6
status = 1;
fprintf(1, ...
'test_IHC3 fails with i = %d, expected_max = %f, blip_max = %f\n', ...
i, expected_maxes(i), blip_maxes(i))
end
if abs(expected_acs(i) - blip_ac(i)) > expected_acs(i)/1e6
if abs(expected_acs(i) - blip_ac(i)) > expected_acs(i)/0.5e6
status = 1;
fprintf(1, ...
'test_IHC3 fails with i = %d, expected_ac = %f, blip_ac = %f\n', ...
Expand Down Expand Up @@ -573,7 +569,7 @@
num_stages = CF.AGC_params.n_stages; % 4
agc_response = zeros(num_stages, CF.n_ch);
for stage = 1:num_stages
agc_response(stage, :) = CF.ears(1).AGC_state.AGC_memory;
agc_response(stage, :) = CF.ears(1).AGC_state(stage).AGC_memory;
end

CF.open_loop = 1; % For measuring impulse response.
Expand Down Expand Up @@ -658,14 +654,14 @@
8000, 3, 8289.883, 3.212
];
case 'do_syn'
results = [ % Preliminary numbers to make test pass...
125, 64, 136.495, 8.648
250, 58, 264.020, 12.909
500, 50, 514.613, 18.320
1000, 40, 1030.546, 27.946
2000, 29, 2038.875, 25.916
4000, 17, 4058.882, 20.959
8000, 2, 8649.710, 4.715
results = [
125, 65, 119.007, 0.238
250, 59, 239.791, 0.942
500, 50, 514.613, 7.249
1000, 40, 1030.546, 30.843
2000, 29, 2038.875, 22.514
4000, 17, 4058.882, 7.691
8000, 4, 7925.624, 1.935
];
end

Expand Down Expand Up @@ -723,7 +719,6 @@
end

%%
% Test: Plot spatial response to match Figure 19.7
if do_plots % Plot final AGC state
figure
plot(agc_response')
Expand Down

0 comments on commit 17f5564

Please sign in to comment.