%simulation for carl hopkins
%
%--generate convolution filters
%--choose input waveform and stim orientation
%--generate receptor outputs -- POISSON output
%--combine to generate small cell outputs

%--This version has full excitatory distribuition!
%--This version has NO head/tail assymetry

%make sure each run is random
rand('state',sum(100*clock));
%figure(4);clf;
figure(2);clf;

%get the excitatory delays in mSec -- 205 entries
fdir = 'C:\Documents and Settings\Bruce Land\My Documents\Matlab\CarlHopkins\HopkinsELA\Hopkins\Feb2006\';
fname = 'ExciteDelay.txt';
load([fdir,fname]);
%====================================================
%define construction parameters
buildEOD = 1;
buildFilters = 1;
buildSmallCells = 1;
%====================================================
%define time
dt = 1e-6; %sec
endTime = 0.004; %sec
t=0:dt:endTime;
%sample frequency
Fs = 1/dt ;

%Initialize Receptor state:
nReceptor = 10;
% vrest = 0;
% rRin = 1 ;
%rVoltage(1:nReceptor) = 0; %resting potential
%rVout(1:nReceptor,length(t))=0;
%rThreshold(1:nReceptor) = 10; %millivolts
%receptor cell time constant
% rTau(1:nReceptor) = 1e-3; %second

%refractory period from Bell and Grant, JNeurosci 1989,p1037
rRefract(1:nReceptor) = 1.2e-3; %second

rIsRefract(1:nReceptor)= 1;
rLastSpike(1:nReceptor) = -inf;
%assign to quadrant 1-4
rPosition = 1 + mod([1:nReceptor],4);

%output spike array initialization
clear rSpike
rSpike= cell(nReceptor,1); %uint8(zeros(length(t),nReceptor));
sVsum = zeros(1,length(t)); %sum of all spikes
%====================================================
%define nReceptor filters, then filter the EOD
%to produce input currents rIin
%Head filters 200 to 4000 Hz
%tail filters 2000+/-500 Hz
if buildFilters
    for i = 1:nReceptor

        %         if rPosition(i)<=2 %the head end
        %             logcf = log(1500) + randn * log(1.75);
        %             cf= exp(logcf);
        %             if cf<300
        %                 cf=1500;
        %             end;
        %             if cf>Fs/2
        %                 cf=1500;
        %             end;
        %         else %the tail end
        %             logcf = log(2000) + randn * log(1.2);
        %             cf= exp(logcf);
        %             if cf<300
        %                 cf=2000;
        %             end;
        %             if cf>Fs/2
        %                 cf=2000;
        %             end;
        %         end
        %MOD for only one set of filters
        logcf = log(2000) + randn * log(1.2);
        cf= exp(logcf);
        if cf<300
            cf=2000;
        end;
        if cf>Fs/2
            cf=2000
        end
        %lower Q for filters
        [b,a]=butter(1,[(cf-cf*0.8)/(Fs/2) (cf+300)/(Fs/2)]);
        %originally 250/300 for low/high width
        bi{i}=b;
        ai{i}=a;
    end
end

%====================================================
%Initialize Small cell state:
if buildSmallCells
    nSmall = 20;
    vrest = 0;
    sRin = 10 ;
    nSmallInputs = 2;
    sThreshold(1:nSmall) = 5; %millivolts
    %cell time constant
    sTau(1:nSmall) = 2e-5; %second
    %refractory period
    sRefract(1:nSmall) = 5e-4; %second
    %Amplitude scaling for input currents
    sAmpInputs = 1 * ones(nSmall,nSmallInputs);
    %time const for input currents to small cells
    sTauInputs(:,1) = 1e-3 * ones(nSmall,1); %second -- inhib
    sTauInputs(:,2:nSmallInputs) = ...
        2e-4 * ones(nSmall,nSmallInputs-1); %second -- excit
    %get random inputs to each small cell
    sInputList = randint(nSmall,nSmallInputs,[1,nReceptor]);
    %gaussian delays
    %inhibitory mean=0.056e-3 std=.011e-3 min=1e-5
    %excitatory mean=0.39e-3 std=.052e-3 min=1e-5
    %inhib
    sInputDelays(:,1) = .002 + ...
        max(0.125e-3*randn(nSmall,1)+0.200e-3, 1e-5) ;
    %excite
    %sInputDelays(:,2:nSmallInputs) = .002 + ...
    %    max(0.052e-3*randn(nSmall,1)+0.386e-3, 1e-5) ;
    sInputDelays(:,2:nSmallInputs) = .002 + ...
        1e-3 * ExciteDelay(randint(nSmall,1,[1,length(ExciteDelay)]));
end

%===Start of per stim code===============================
%define stimulus orientations
%order frontleft, frontright, backleft, backright
Stimulus(1,:) = [1 1 -1 -1] ;%front/back pos front
Stimulus(2,:) = [-1 -1 1 1] ;%front/back neg front
Stimulus(3,:) = [1 -1 1 -1] ;%left/right pos left
Stimulus(4,:) = [-1 1 -1 1] ;%left/right neg left
Stimulus(5,:) = [ 1 1 1 1 ] ; %pos outside
Stimulus(6,:) = [ -1 -1 -1 -1 ] ; %pos outside
%set stim type to 1-4 (as above)
StimType = 1;
EODampList = [.002 .005 .01 .02, .05 .08 .1 .15 .2 .3 .5 1 2 5];
%====================================================
%make an EOD
% EOD1amp = 1;
% EOD1dur1 = .0005;
% EOD1dur2 = .001;
% EOD2amp = 0.5;
% %concantenating two 1/2 sin cycles
% oneT = 0:dt:(EOD1dur1+EOD1dur2);
% EOD = EOD1amp*sin(2*pi/(2*EOD1dur1)*oneT).*(oneT<=EOD1dur1) + ...
%     EOD2amp*sin(2*pi/(2*EOD1dur2)*(oneT+EOD1dur1))...
%     .*(oneT>EOD1dur1) ;
% %pad with zeros at the beginning and end
% %so it is the same length as the vector t
% EOD = [zeros(1,100),...
%     EOD,...
%     zeros(1,length(t)-length(EOD)-100)];
if buildEOD
    %make a square wave
    sqDur = 2e-4;%duration in mSec
    EOD = zeros(1,length(t));
    EOD(100:fix(100+sqDur/dt)) = 1; %amplitude 0.1 is small 1 is big
end

%====================================================
%now actually filter
for i = 1:nReceptor
    rIin(:,i) =  filter(bi{i},ai{i},EOD);
end
%apply the stimulus type the the filtered current
rIin =  rIin .* ...
    repmat(Stimulus(StimType,rPosition),length(t),1);
%apply a perfect rectifier
rIin = rIin .* (rIin>0); %rectify

%====================================================
%Compute receptor outputs accurate to 1 microsec
%update voltage
tindex=1;
for time=t
    %update the voltage
    %     rVoltage = rVoltage + ...
    %         (dt./rTau) .* ...
    %         (rRin .* rIin(tindex,:) - (rVoltage - vrest));

    %generate a spike
    %note that rIsRefract is zero after a spike
    tspike = rand(1,nReceptor) < .05*(rIin(tindex,:).*rIsRefract);
    %record it for later
    %rSpike(tindex,:)=uint8(tspike);
    for i=1:nReceptor
        if tspike(i)==1
            rSpike{i} = [ [rSpike{i}],time] ;
        end
    end
    %record the last spike time for refractory period
    rLastSpike(tspike)=time;
    %reset the spiking cells
    %rVoltage(tspike) = vrest;
    %enforce refractory period
    rIsRefract = (rLastSpike+rRefract) < time;

    %sum up the APs as if they are sq waves
    sVsum(tindex) = sVsum(tindex) + sum(tspike);
    %update output array index
    %rVout(:,tindex) = rVoltage;
    tindex = tindex + 1;
end

figure(1);clf;
col =  [1 0 0
    0 1 0
    0 0 1
    0 1 1];

subplot(3,1,1)
for i=1:nReceptor
    plot(t,rIin(:,i),'color',col(rPosition(i),:))
    hold on
end
title('input currents')
plot(t,EOD,'k')
% subplot(3,1,2)
% for i=1:nReceptor
%     plot(t,rVout(i,:),'color',col(rPosition(i),:))
%     hold on
% end
%plot(t,EOD,'k')
title('Receptor voltages')

subplot(3,1,2)
%make a dot-plot of all the receptors
for i=1:nReceptor
    plot(rSpike{i},i*ones(1,length(rSpike{i})),...
        '.','color',col(rPosition(i),:));
    hold on
end
title('Receptors r=FL; g=FR, b=BL, y=BR')
set(gca,'xlim',[0,endTime])

subplot(3,1,3)
gsd = 10e-5; %resolution of averaging window
tg = -3*gsd:dt:3*gsd;
gg = exp(-(tg.^2/gsd^2));
ss = conv(sVsum,gg);
ssr = ss(fix(length(tg)/2):end-fix(length(tg)/2)-1);
% [b,a]=butter(4,.01);
% ss = filtfilt(b,a,sVsum);
%
plot(t,ssr,'k')
%bar(sum(reshape(sVsum(1:4000),20,200)))

%====================================================
%Compute small cell outputs accurate to 1 microsec
%init output arrays
clear sSpike
sSpike = cell(nSmall,1);;
sVoltage(1:nSmall) = 0;
sVout(1:nSmall,length(t)) = 0; %all voltages
sVsum = zeros(1,length(t)); %sum of all spikes
sVpart = zeros(4,length(t)); %sum of spikes with each origin

for i=1:nSmall
    current = zeros(1,length(t));
    for j=1:nSmallInputs
        %define alpha function: sAmpInputs*exp(-t/sTauInputs)
        %carry it out for 4 time constants
        formTime = 0:dt:4*sTauInputs(i,j);
        form = sAmpInputs(i,j)*exp(-formTime/sTauInputs(i,j));
        if j==1
            form = -10*form; %inhibitory input
        end
        %build input spikes (sReceptor shifted by delays)
        %add the random delay times to the receptor
        %spike output train chosen by sImputList
        delayedSpikes = ...
            rSpike{sInputList(i,j)} + sInputDelays(i,j);
        %add a copy of the alpha function into the current
        %for every spike
        for k=1:length(delayedSpikes)
            %make this longer than it has to be
            delayedform = ...
                [zeros(1,fix(delayedSpikes(k)/dt)), form, zeros(1,length(t))];
            %then truncate
            delayedform = delayedform(1:length(t));
            current = current + delayedform ;
        end
    end %loop thru all inputs to one small cell
    %at this point we have the input current to one Smallcell,
    %so now we can generate Small cell voltages
    %update voltage
    sLastSpike = -inf;
    tindex=1;
    for time=t
        %update the voltage
        sVoltage(i) = sVoltage(i) + ...
            (dt./sTau(i)) .* ...
            (sRin .* current(tindex) - (sVoltage(i) - vrest));

        %generate a spike -- add threshold noise?
        tspike = (sVoltage(i) >= sThreshold(i));
        %record it for later
        %rSpike(tindex,:)=uint8(tspike);

        if tspike==1
            sSpike{i} = [ [sSpike{i}],time] ;
            %record the last spike time for refractory period
            sLastSpike=time;
            %reset the spiking cells
            sVoltage(i) = vrest;
        end
        %enforce refractory period
        sVoltage(i) = sVoltage(i) .*...
            ((sLastSpike + sRefract(i)) < time);

        %sum up the APs as if they are sq waves
        sVsum(tindex) = sVsum(tindex) + tspike;
        %sum up each origin
        sVpart(rPosition(sInputList(i,2)),tindex)= ...
            sVpart(rPosition(sInputList(i,2)),tindex) + tspike;
        %update output array index
        sVout(i,tindex) = sVoltage(i);

        tindex = tindex + 1;
    end %time loop
    figure(2)
    subplot(4,1,1)
    plot(t,current,'color',col(rPosition(sInputList(i,2)),:))
    title('Small cell input currents')
    hold on
    subplot(4,1,2)
    plot(t,sVout(i,:),'color',col(rPosition(sInputList(i,2)),:))
    title('Small cell intracellular voltages')
    hold on
end %loop thru all small cells
plot(t,EOD,'k')
%====================================================
subplot(4,1,3)
%make a dot-plot of all the receptors
for i=1:nSmall
    plot(sSpike{i},i*ones(1,length(sSpike{i})),...
        '.','color',col(rPosition(sInputList(i,2)),:));
    hold on
end
title('Small Cells -- color of excitatory receptor input')
set(gca,'xlim',[0,endTime])
%====================================================
%sum up all the small cell APs
subplot(4,1,4)
gsd = 10e-5; %resolution of averaging window
tg = -3*gsd:dt:3*gsd;
gg = exp(-(tg.^2/gsd^2));
ss = conv(gg,sVsum);
sss = ss(fix(length(tg)/2):end-fix(length(tg)/2)-1);
%  [b,a]=butter(4,.010);
%  ss = filtfilt(b,a,sVsum);
plot(t,sss,'k')
hold on
for j=1:4
    ss = conv(gg,sVpart(j,:));
    ss = ss(fix(length(tg)/2):end-fix(length(tg)/2)-1);
    %ss = filtfilt(b,a,sVpart(j,:));
    plot(t,ss...
        ,'color',col(j,:));
end
title('Summed Small Cell APs')
%====================================================
%format the total so that it looks like the physiology
figure(3);clf;
plot(t,-(sss)) %+ssr
