%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
%--This version plots MULTIPLE amp/resonse curves
clear all
figure(2);clf;
figure(1);clf;
%figure(3);clf;


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

%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;
nFish = 1;
OneFish = 1; %Set to zero for multiple fish
%====================================================
%define time
dt = 1e-6; %sec
endTime = 0.002; %sec originally .004
t=0:dt:endTime;
%sample frequency
Fs = 1/dt ;

%Initialize Receptor state:
nReceptor = 100;
% 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

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

for FishNum=1:nFish
    if FishNum>1 & OneFish==1
        buildFilters = 0;
        buildSmallCells = 0;
    end
    %====================================================
    %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 = 800;
        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
        %==inhib==
        %inhibitory mean=0.056e-3 std=.011e-3 from Friedman etal
        inhibMean = .2e-3 ;%+ .15e-3 ; %0.056e-3 add up to .3 mSec for syn delay
        inhibStd  = .1e-3;  %  0.011e-3;% + 0.03e-3 ;
        propDelay = 0; %set to .002 for better agreement with real data
        sInputDelays(:,1) = propDelay + ...
            max(inhibStd*randn(nSmall,1)+inhibMean, 1e-5) ;
        %==excite==
        %excitatory mean=0.39e-3 std=.052e-3 min=1e-5
        %sInputDelays(:,2:nSmallInputs) = .002 + ...
        %    max(0.052e-3*randn(nSmall,1)+0.386e-3, 1e-5) ;
        %drawn from the table of measured values
        sInputDelays(:,2:nSmallInputs) = propDelay + ...
            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 = 5;
    %EODampList = [.002 .005 .01 .02, .03 .05 .08 .1 .2 .3 .5 1 2 5];
    StimTypeList = [1 1 1 1 1  2 2 2 2 2  3 3 3 3 3  4 4 4 4 4];
    %StimTypeList = [1 1 1 1 1 1    3 3 3 3 3 3];
    
    EODampList = [.01 .05 .2 .5 1  .01 .05 .2 .5 1  .01 .05 .2 .5 1  .01 .05 .2 .5 1  ];
    sqDur = 3e-4;%duration in mSec
    EODdurList = [sqDur sqDur sqDur sqDur sqDur...
        sqDur sqDur sqDur sqDur sqDur ...
        sqDur sqDur sqDur sqDur sqDur ...
        sqDur sqDur sqDur sqDur sqDur ];

%     EODampList = [.1 .05 .5 .05 .5 .1 .05 .5 .05 .5 .1 .05 .5 .05 .5 .1 .05 .5 .05 .5 ];
%     EODdurList = [.0003 .0002 .0002 .001 .001 ...
%         .0003 .0002 .0002 .001 .001 ...
%         .0003 .0002 .0002 .001 .001 ...
%         .0003 .0002 .0002 .001 .001];

%     EODampList = 0.1 * [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
%     EODdurList = [.0001 .0002 .0003 .0005 .001 ...
%         .0001 .0002 .0003 .0005 .001 ...
%         .0001 .0002 .0003 .0005 .001 ...
%         .0001 .0002 .0003 .0005 .001];

    %binary occurance of a given cell firing
    sSpike = zeros(nSmall,length(EODampList)) ;

    for EODindex=1:length(EODampList)
        EODindex
        %init refactory period
        rIsRefract(1:nReceptor)= 1;
        rLastSpike(1:nReceptor) = -inf;

        %output spike array initialization
        clear rSpike
        rSpike= cell(nReceptor,1); %uint8(zeros(length(t),nReceptor));
        sVsum = zeros(1,length(t)); %sum of all spikes
        %====================================================
        %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

            EOD = zeros(1,length(t));
            EOD(100:fix(100+EODdurList(EODindex)/dt)) = EODampList(EODindex); %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(StimTypeList(EODindex),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

        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,EODindex) = 1 ;
                    %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')
        %====================================================
        %====================================================
        %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);
        %figure(2)
        %col = jet(length(EODampList));
        %plot(sss,'color',col(EODindex,:))
        %drawnow
        %hold on
        %find the max and save it
        maxSmallResponse(EODindex) = max(sss);
        %  [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')
        %====================================================
        %         figure(2)
        %         col = jet(length(EODampList));
        %         plot(sSpike, 'color', col(EODindex,:) )
        %         set(gca,'ylim',[0 1.2])
        %         hold on
        %         drawnow
    end %EOD amp list

    figure(2)
    for i=1:length(EODampList)
        plot(sSpike(:,i)*.8+i)
        hold on
    end


    figure(1)
    %col = jet(nFish);
    %col = min(gray(nFish)+0.1, 1);
    group1 = find(StimTypeList==1);
    group2 = find(StimTypeList==2);
    group3 = find(StimTypeList==3);
    group4 = find(StimTypeList==4);
    semilogx(EODampList(group1),maxSmallResponse(group1),'-or')
    hold on
    semilogx(EODampList(group2),maxSmallResponse(group2),'-og')
    semilogx(EODampList(group3),maxSmallResponse(group3),'-ob')
    semilogx(EODampList(group4),maxSmallResponse(group4),'-oc')

    title([  ...
        '  inhibM=',num2str(inhibMean),'  inhibSt=',num2str(inhibStd),...
        '  N=',num2str(nFish),'  OneFish=',num2str(OneFish), ...
        '  nSmall=',num2str(nSmall)])

    text(.01,max(maxSmallResponse)-5,num2str(StimTypeList,'%5d '),'FontSize',8)
    text(.01,max(maxSmallResponse)-10,num2str(1000*EODdurList,'%5.1f'),'FontSize',8)
    text(.01,max(maxSmallResponse)-15,num2str(EODampList,'%5.2f'),'FontSize',8)
    drawnow
    hold on

    %similarity matrix
    %figure(3)
    for ii=1:length(EODampList)
        for jj=1:length(EODampList)
            dist(ii,jj) = sum(abs(sSpike(:,ii)-sSpike(:,jj)));
        end
    end
    dist
    
    % ----------------------------------------------------------------------------------
    %    MULTIDIMENSIONAL SCALING in matlab by Mark Steyvers 1999
    %      needs optimization toolbox
    %
    %Modified by Bruce Land
    %--Data via globals to anaylsis programs
    %--3D plotting with color coded groups
    %--Mapping of MDS space to spike train 'temporal profiles' as described in
    %Aronov, et.al. "Neural coding of spatial phase in V1 of the Macaque" in
    %press J. Neurophysiology
    % -------------------------------------------------------------------------
    % ---------
    %no 			= 10;						% size of input matrix
    %A  			= rand( no , no );	% dissimilarities
    %load 'C:\Documents and Settings\bruce land\My Documents\Matlab\Damian\spdtable06dec';
    %load 'C:\Documents and Settings\bruce land\My Documents\Matlab\Damian\spctable06dec';

    % no = 3;
    % A = [ 0 3 4 ; ...
    %       3 0 5; ...
    %       4 5 0 ];
    no = length(EODampList);

    ngrp = 4;
    grp{1} = 1:5  ;
    grp{2} = 6:10 ;
    grp{3} = 11:15 ;
    grp{4} = 16:20 ;


    doHull=0; %if 1 does convex hull plot

    A = dist ;
    %anovadim = 1;
    plotdim1 	= 1;		% which dimension on x-axis
    plotdim2 	= 2;		% which dimension on y-axis
    plotdim3    = 3;
    ndims    	=3;		% number of dimension in MDS solution
    R			= 2.0;	% R values in minkowski metric
    maxiter		= 50;		% maximum number of iterations
    conv     	= 0.001;	% convergence criterion
    R1       	= 1;     % 1=Torgeson Young scaling for initial configuration 0 = intial random configuration
    seed     	= 1;		% seed for random number generator
    minoption 	= 1;  	% 1 = minimize stress1    2 = minimize stress2

    % provide some text labels for the stimulus points here
    for i=1:no
        labels{ i } = sprintf( '%d' , i );
    end

    userinput{2} = R;
    userinput{3} = ndims;
    userinput{4} = maxiter;
    userinput{5} = conv;
    userinput{6} = 1;         % 0=no 1=yes, printed comments
    userinput{7} = 0;
    userinput{8} = 0;
    userinput{9} = R1;
    userinput{10}= 0;
    userinput{11}= seed;
    userinput{13}= 1;         % 0=do not symmetrize input matrix % 1=do symmetrize
    userinput{14}= 0;
    userinput{15}= minoption;
    
    % ---------------------------------------------------------------------------------
    %    CALL THE MDS ROUTINE
    [ Config,DHS,DS,DeltaS,Stress1,StressT1,Stress2,StressT2,Rs,RsT] = ...
        mds( userinput,A );
    % ---------------------------------------------------------------------------------

    % ---------------------------------------------------------------------------------
    %                            Create the Shepard Plot
    % ---------------------------------------------------------------------------------
    figure( 3 );clf;
    plot( DS , DeltaS , 'r*' , DHS , DeltaS , '-gs' );
    grid on;
    axis square;
    axis tight;
    ylabel( 'dissimilarities' );
    xlabel( 'distances' );
    title( sprintf( 'stress1=%1.4f stress2=%1.4f Rs=%1.4f (N=%d)' , Stress1 , Stress2 , Rs , no ) , 'FontSize' , 8 );

    % ---------------------------------------------------------------------------------
    %                            Create the plot with stimulus coordinates
    % ---------------------------------------------------------------------------------
    figure( 4 );clf;

    if (ndims==1)

        plot( Config(:,1) , Config(:,1) , 'r*' );
        grid on;
        axis square;
        xlabel( 'dimension 1' );
        ylabel( 'dimension 1' );

        for i=1:no
            text( Config(i,1)+0.1 , Config(i,1) , labels{i}  );
        end;
    elseif (ndims>=2)
        %plot( Config(:,plotdim1) , Config(:,plotdim2) , 'r*' );
        map = [1 0 0; 0 1 0; 0 0 1; 0 1 1; 0 0 0];
        for i=1:ngrp
            pts = grp{i};
            col = map(i,:);
            plot( Config(pts,plotdim1),Config(pts,plotdim2) ,...
                'd','color',col,'markerfacecolor',col,'markersize',10 );
            hold on
            if i<5
                k = convhull(Config(pts,plotdim1),Config(pts,plotdim2));
                plot(Config(k+pts(1)-1,plotdim1),Config(k+pts(1)-1,plotdim2),'color',col)
            end
        end

        grid on;
        axis equal;
        xlabel( sprintf( 'dimension %d' , plotdim1 ));
        ylabel( sprintf( 'dimension %d' , plotdim2 ));
        %
         title([  ...
        '  inhibM=',num2str(inhibMean),'  inhibSt=',num2str(inhibStd),...
        '  N=',num2str(nFish),'  OneFish=',num2str(OneFish), ...
        '  nSmall=',num2str(nSmall)])

        text(min(Config(1:length(EODampList),plotdim1)),...
            max(Config(1:length(EODampList),plotdim2)),...
            num2str(StimTypeList,'%6d '),'FontSize',8)
        text(min(Config(1:length(EODampList),plotdim1)),...
            max(Config(1:length(EODampList),plotdim2))-.15,...
            num2str(1000*EODdurList,'%6.1f'),'FontSize',8)
        text(min(Config(1:length(EODampList),plotdim1)),...
            max(Config(1:length(EODampList),plotdim2))-.30,...
            num2str(EODampList,'%5.2f'),'FontSize',8)
    %
        for i=1:no
            text( Config(i,plotdim1)+0.1 , Config(i,plotdim2) , labels{i} );
        end;
    end

    if (ndims==3)
        figure(5);clf;

        az = 60;
        %subplot(1,2,1)
        map = [1 0 0; 0 1 0; 0 0 1;0 1 1; 0 0 0];
        for i=1:ngrp
            pts = grp{i};
            col = map(i,:);
            plot3( Config(pts,plotdim1),Config(pts,plotdim2),Config(pts,plotdim3) ,...
                '*','color',col );
            hold on
        end
        min3 = get(gca,'zlim');
        for i=1:ngrp
            pts = grp{i};
            col = map(i,:);
            line([Config(pts,plotdim1) Config(pts,plotdim1)]',...
                [Config(pts,plotdim2) Config(pts,plotdim2)]',...
                [Config(pts,plotdim3) min3(1)*ones(size(Config(pts,plotdim3)))]',...
                'color',col,'linewidth',2  )
        end

        grid on;
        box on;
        axis equal;
        xlabel( sprintf( 'dimension %d' , plotdim1 ));
        ylabel( sprintf( 'dimension %d' , plotdim2 ));
        zlabel( sprintf( 'dimension %d' , plotdim3 ));

        for i=1:no
            text( Config(i,plotdim1)+0.1 ,Config(i,plotdim2),Config(i,plotdim3),...
                labels{i},'fontname','helvetica','fontweight','bold' );
        end;
        view(az,15)
        set(gca,'projection','perspective')
        rotate3d on
    end
end