%% License Info
%A Model of Feedforward, Global, and Lateral Inhibition in the Locust Visual System.

%This model examines the architecture and function of inhibitory mechanisms in the 
%visual system of locusts, namely those involved in the processing of inputs to a
%key looming-sensitive neuron, the lobula giant movement detector (LGMD).

%Copyright (c) 2026, Erik Olson, Travis Wiens, and Jack Gray

%CITATION:
%When using the model code for scientific publications, cite the following work:
%Olson EGN, Wiens TK, Gray JR. A model of feedforward, global, and lateral inhibition 
%in the locust visual system predicts responses to looming stimuli. Biol Cybern. 2021 
%Jun;115(3):245-265. doi: 10.1007/s00422-021-00876-8. Epub 2021 May 16. PMID: 33997912.

%This program is free software: you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation, either version 3 of the License, or
%(at your option) any later version.

%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.

%You should have received a copy of the GNU General Public License
%along with this program.  If not, see <https://www.gnu.org/licenses/>.

%Contact: erik.olson@usask.ca

%% Code
function [f_output,f_all,theta_q,raw,t_raw,V_trace] = FFITrialFunction(parFile,smoothing,l_over_v,inhibition_mod)
%Testing FFI output in response to looming stimuli centered on an
%ommatidium, across all possible offsets with respect to the symmetry of 
%the FFI neuron distribution pattern.  Inputs:
%parfile = name of the parameter file containing model parameters
%smoothing = width of the Gaussian window [ms] used to generate firing rate
%values from the spike timing histogram
%l_over_v = stimulus half-height to velocity ratio [s] used to determine
%stimulus size
%inhibition_mod = vector of two values which can override lateral and
%global inhibitory gains, respectively.  To avoid override, set to NaN.
%Outputs:
%f_output = population firing rate [spikes/s] summed across the modelled
%FFI neurons, averaged across trials with all possible offsets.
%f_all = population firing rate [spikes/s]  summed across the modelled FFI 
%neurons, with each column corresponding to a unique offset.
%theta_q = visual subtense angles [deg] corresponding to firing rate
%values in f_output and f_all.
%raw = raw spiketrain data for all trials
%t_raw = corresponding time values [s] for "raw" and "V_trace"
%V_trace = voltage traces from modelled FFI neurons for a single looming
%trial, centered on the eye

%keep only the last 1s of data from each loom
cut = 1;

%try
[f,t_f,raw,t_raw,V_trace] = subIterationFFI(parFile,smoothing,l_over_v,0,0,inhibition_mod,cut);

results = NaN(size(f,1),19);
raw_results = NaN(size(raw,1),19);
results(:,1) = f;
raw_results(:,1) = raw;

h_offsets = [0,1,1,2,2,2,3,3,3,3,4,4,4,4,4,5,5,5,6];
v_offsets = [0,0,1,0,1,2,0,1,2,3,0,1,2,3,4,0,1,2,0];

parfor i = 2:19
    h_offset = h_offsets(i);
    v_offset = v_offsets(i);
    [f,~,raw] = subIterationFFI(parFile,smoothing,l_over_v,h_offset,v_offset,inhibition_mod,cut);
    results(:,i) = f;
    raw_results(:,i) = raw;
end

load('Wang_et_al_2018_FFI_Looming_Data.mat','theta_exp_20ms');
theta_exp = theta_exp_20ms;

f = results;
raw = raw_results;

load(parFile,'v','d');
l = v*l_over_v;
d = cut*v;
theta = (2*atan(l./(d-t_f*v)))*180/pi;

theta_q = linspace(theta(1),theta_exp(end),1001)';

f2 = interp1(theta,f,theta_q);

f_output = mean(f2,2);
f_all = f2;

fprintf('Trial complete with l/v = %.4f \n',l_over_v);
end