Move SYN constants from code into default parameters

Moving a bunch of magic numbers out of the SYN_Design code and into the default SYN_params, so they are modifiable by user code.  And various other code tweaks and naming improvements (e.g. changed z to r for release rate; short variable names to facilitate labeling on doc diagrams).

matlab/CARFAC_Design.m

@@ -150,13 +150,19 @@
 
 if num_args < 6  % Default the SYN_params.
   n_classes = 3;  % Default.  Modify params and redesign to change.
-  % Parameters could generally have columns if channel-dependent.
+  % Parameters could generally have columns if class-dependent.
   CF_SYN_params = struct( ...
     'do_syn', do_syn, ...  % This may just turn it off completely.
     'fs', fs, ...
     'n_classes', n_classes, ...
+    'n_fibers', [50, 35, 25], ...
+    'spont_rates', [50, 6, 1], ...
+    'sat_rates', 200, ...
+    'sat_reservoir', 0.2, ...
+    'v_width', 0.02, ...
     'tau_lpf', 0.000080, ...
-    'reservoir_tau', 0.020);
+    'reservoir_tau', 0.020, ...
+    'agc_weights', [1.2, 1.2, 1.2]');  % column
 end
 
 % first figure out how many filter stages (PZFC/CARFAC channels):
@@ -577,55 +583,37 @@
 end
 
 n_ch = length(pole_freqs);
-n_cl = SYN_params.n_classes;
+n_classes = SYN_params.n_classes;
 
-sat_rates = 200;
-switch n_cl  % Just enough to test that both 2 and 3 can work.
-  case 2
-    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
-    spont_rates = [50, 6, 1];
-    n_fibers = [50, 35, 25];
-    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);
+v_widths = SYN_params.v_width * ones(1, n_classes);
 
 % 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;
+p1 = SYN_params.sat_rates / fs;
+p0 = SYN_params.spont_rates / fs;
 
+w1 = SYN_params.sat_reservoir;
 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;
+% Assume the sigmoid is switching between 0 and 1 at 50% duty cycle, so
+% normalized mean value is 0.5 at saturation.
+s1 = 0.5;
+r1 = s1*w1;
+% solve q1 = a1*r1 for gain coeff a1:
+a1 = q1 ./ r1;
+% solve p1 = a2*r1 for gain coeff a2:
+a2 = p1 ./ r1;
 
 % Now work out how to get the desired spont.
-z0 = p0 ./ a2;
-q0 = z0 .* a1;
+r0 = p0 ./ a2;
+q0 = r0 .* a1;
 w0 = 1 - q0;
-s0 = z0 ./ w0;
+s0 = r0 ./ w0;
 % Solve for (negative) sigmoid midpoint offset that gets s0 right.
 offsets = log((1 - s0)./s0);
 
@@ -634,14 +622,15 @@
 % Copy stuff needed at run time into coeffs.
 SYN_coeffs = struct( ...
   'n_ch', n_ch, ...
-  'n_classes', n_cl, ...
-  'n_fibers', ones(n_ch,1) * n_fibers, ...  % Synaptopathy comes in here.
+  'n_classes', n_classes, ...
+  'n_fibers', ones(n_ch,1) * SYN_params.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.
-  'spont_p', spont_p, ...
+  'agc_weights', SYN_params.agc_weights, ... % For making a nap out to agc in.
+  'spont_p', spont_p, ... % used only to init the output LPF
+  'spont_sub', (SYN_params.n_fibers .* spont_p) * SYN_params.agc_weights, ...
   'res_lpf_inits', q0, ...
   'res_coeff', 1 - exp(-1/(SYN_params.reservoir_tau * fs)), ...
   'lpf_coeff', 1 - exp(-1/(SYN_params.tau_lpf * fs)));

matlab/CARFAC_SYN_Step.m

@@ -6,21 +6,26 @@
 % 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.
-s = 1 ./ (1 + exp(-x));
-release_rates = (1 - state.reservoirs) .* s;  % z = w*s
+s = 1 ./ (1 + exp(-x));  % Between 0 and 1; positive at rest.
+q = state.reservoirs;  % aka 1 - w, between 0 and 1; positive at rest.
+r = (1 - q) .* s;  % aka w*s, between 0 and 1, proportional to release rate.
 
-% Smooth once with LPF (receptor potential was already smooth):
+% Smooth once with LPF (receptor potential was already smooth), after
+% applying the gain coeff a2 to convert to firing prob per sample.
 state.lpf_state = state.lpf_state + coeffs.lpf_coeff * ...
-  (coeffs.a2 .* release_rates - state.lpf_state);  % this is firing probs.
+  (coeffs.a2 .* r - 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 .* firing_probs;
 
 % Feedback, to update reservoir state q for next time.
-state.reservoirs = state.reservoirs + coeffs.res_coeff .* ...
-  (coeffs.a1 .* release_rates - state.reservoirs);
+state.reservoirs = q + coeffs.res_coeff .* (coeffs.a1 .* r - q);
 
-% Make an output that resembles ihc_out, to go to agc_in.
-syn_out = (coeffs.n_fibers .* (firing_probs - coeffs.spont_p)) * coeffs.agc_weights;
+% Make an output that resembles ihc_out, to go to agc_in (collapse over classes).
+% Includes synaptopathy's presumed effect of reducing feedback via n_fibers.
+% But it's relative to the healthy nominal spont, so could potentially go
+% a bit negative in quiet is there was loss of high-spont or medium-spont units.
+% The weight multiplication is an inner product, reducing n_classes
+% columns to 1 column.
+syn_out = (coeffs.n_fibers .* firing_probs) * coeffs.agc_weights - coeffs.spont_sub;