From 17f55640ea7504f6cadbb8f0b34b5babdfcb6e05 Mon Sep 17 00:00:00 2001 From: Richard Lyon Date: Thu, 19 Sep 2024 11:33:03 -0700 Subject: [PATCH] Reparameterized, better design phase 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. --- matlab/CARFAC_Design.m | 81 ++++++++++++++++++++++++++-------------- matlab/CARFAC_IHC_Step.m | 7 ++-- matlab/CARFAC_Init.m | 5 +-- matlab/CARFAC_SYN_Step.m | 27 +++++++------- matlab/CARFAC_Test.m | 65 +++++++++++++++----------------- 5 files changed, 104 insertions(+), 81 deletions(-) diff --git a/matlab/CARFAC_Design.m b/matlab/CARFAC_Design.m index 94fb3d9..2a38187 100644 --- a/matlab/CARFAC_Design.m +++ b/matlab/CARFAC_Design.m @@ -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): @@ -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))); diff --git a/matlab/CARFAC_IHC_Step.m b/matlab/CARFAC_IHC_Step.m index 0c9e133..8af937f 100644 --- a/matlab/CARFAC_IHC_Step.m +++ b/matlab/CARFAC_IHC_Step.m @@ -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 @@ -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 @@ -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. @@ -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 diff --git a/matlab/CARFAC_Init.m b/matlab/CARFAC_Init.m index 889330c..3be9f7e 100644 --- a/matlab/CARFAC_Init.m +++ b/matlab/CARFAC_Init.m @@ -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); diff --git a/matlab/CARFAC_SYN_Step.m b/matlab/CARFAC_SYN_Step.m index a15af03..8ea2b82 100644 --- a/matlab/CARFAC_SYN_Step.m +++ b/matlab/CARFAC_SYN_Step.m @@ -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; diff --git a/matlab/CARFAC_Test.m b/matlab/CARFAC_Test.m index b170dd4..77b782a 100644 --- a/matlab/CARFAC_Test.m +++ b/matlab/CARFAC_Test.m @@ -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, :); @@ -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 @@ -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 @@ -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', ... @@ -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', ... @@ -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. @@ -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 @@ -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')