comparison utilities/UTIL_showMAP.m @ 38:c2204b18f4a2 tip

End nov big change

-function UTIL_showMAP (showMapOptions, paramChanges)
+function UTIL_showMAP (showMapOptions)
 % UTIL_showMAP produces summaries of the output from MAP's mmost recent run
-%  All MAP outputs are stored in global variables and UTIL_showMAP
-%  simply assumes that they are in place.
+%  All MAP_14 outputs are stored in global variables and UTIL_showMAP
+%   simply assumes that they are in place. It uses whatever it there.
 %
-% showMapOptions
+% showMapOptions:
 % showMapOptions.printModelParameters=1; % print model parameters
 % showMapOptions.showModelOutput=1;      % plot all stages output
 % showMapOptions.printFiringRates=1;     % mean activity at all stages
 % showMapOptions.showACF=1;              % SACF (probabilities only)
 % showMapOptions.showEfferent=1;         % plot of efferent activity
-% showMapOptions.surfProbability=0;      % HSR (probability) surf plot
-% showMapOptions.fileName=[];            % parameter filename for plot title
-
+% showMapOptions.surfAN=0;               % AN output surf plot
+% showMapOptions.fileName='';            % parameter filename for plot title
+% showMapOptions.PSTHbinwidth=0.001      % binwidth for surface plots
+% showMapOptions.colorbar=1;             % request color bar if appropriate
+% showMapOptions.view=[0 90];
 % dbstop if warning
 
-global dt ANdt  savedBFlist saveAN_spikesOrProbability saveMAPparamsName...
-    savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...
-    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
-    IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas  ...
-    CNtauGk CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates ...
-    MOCattenuation
-global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams
-global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams
+% Discover results left behind by MAP1_14
+global  savedParamChanges savedBFlist saveAN_spikesOrProbability ...
+    saveMAPparamsName savedInputSignal dt dtSpikes ...
+    OMEextEarPressure TMoutput OMEoutput DRNLoutput...
+    IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
+    IHCoutput ANprobRateOutput ANoutput savePavailable saveNavailable ...
+    ANtauCas CNtauGk CNoutput ICoutput ICmembraneOutput ICfiberTypeRates...
+    MOCattenuation ARattenuation 
+
+% Find parameters created in MAPparams<name>
+global inputStimulusParams OMEParams DRNLParams IHC_cilia_RPParams
+global IHCpreSynapseParams  AN_IHCsynapseParams
+global MacGregorParams MacGregorMultiParams  filteredSACFParams
 global ICrate
 
 
 restorePath=path;
 addpath ( ['..' filesep 'utilities'], ['..' filesep 'parameterStore'])
 
-if nargin<1
-    showMapOptions=[];
-end
+if nargin<1,  showMapOptions=[]; end
+paramChanges=savedParamChanges; 
+
 % defaults (plot staged outputs and print rates only) if options omitted
 if ~isfield(showMapOptions,'printModelParameters')
     showMapOptions.printModelParameters=0; end
-if ~isfield(showMapOptions,'showModelOutput'),showMapOptions.showModelOutput=1;end
-if ~isfield(showMapOptions,'printFiringRates'),showMapOptions.printFiringRates=1;end
+if ~isfield(showMapOptions,'showModelOutput'),showMapOptions.showModelOutput=0;end
+if ~isfield(showMapOptions,'printFiringRates'),showMapOptions.printFiringRates=0;end
+if ~isfield(showMapOptions,'surfAN'),showMapOptions.surfAN=0;end
+if ~isfield(showMapOptions,'ICrates'),showMapOptions.ICrates=0;end
 if ~isfield(showMapOptions,'showACF'),showMapOptions.showACF=0;end
 if ~isfield(showMapOptions,'showEfferent'),showMapOptions.showEfferent=0;end
-if ~isfield(showMapOptions,'surfProbability'),showMapOptions.surfProbability=0;end
+if ~isfield(showMapOptions,'PSTHbinwidth'),showMapOptions.PSTHbinwidth=0.001;end
 if ~isfield(showMapOptions,'fileName'),showMapOptions.fileName=[];end
-if ~isfield(showMapOptions,'surfSpikes'),showMapOptions.surfSpikes=0;end
-if ~isfield(showMapOptions,'ICrates'),showMapOptions.ICrates=0;end
+if ~isfield(showMapOptions,'colorbar'),showMapOptions.colorbar=1;end
+if ~isfield(showMapOptions,'view'),showMapOptions.view=[-25 80];end
+if ~isfield(showMapOptions,'ICrasterPlot'),showMapOptions.ICrasterPlot=0;end
 
 %% send all model parameters to command window
 if showMapOptions.printModelParameters
     % Read parameters from MAPparams<***> file in 'parameterStore' folder
     %  and print out all parameters
-    if nargin>1
         cmd=['MAPparams' saveMAPparamsName '(-1, 1/dt, 1, paramChanges);'];
         eval(cmd);
-    else
-        cmd=['MAPparams' saveMAPparamsName '(-1, 1/dt, 1);'];
-        eval(cmd);
-        disp(' no parameter changes found')
-    end
-end
+end
+
+% ignore zero elements in signal
+signalRMS=mean(savedInputSignal(find(savedInputSignal)).^2)^0.5;
+signalRMSdb=20*log10(signalRMS/20e-6);
+nANfiberTypes=length(ANtauCas);
 
 %% summarise firing rates in command window
 if showMapOptions.printFiringRates
     %% print summary firing rates
     fprintf('\n')
     disp('summary')
     disp(['AR: ' num2str(min(ARattenuation))])
     disp(['MOC: ' num2str(min(min(MOCattenuation)))])
-    nANfiberTypes=length(ANtauCas);
     if strcmp(saveAN_spikesOrProbability, 'spikes')
         nANfibers=size(ANoutput,1);
         nHSRfibers=nANfibers/nANfiberTypes;
         duration=size(TMoutput,2)*dt;
         disp(['AN(HSR): ' num2str(sum(sum(ANoutput(end-nHSRfibers+1:end,:)))/...

 
 
 %% figure (99): display output from all model stages
 if showMapOptions.showModelOutput
     plotInstructions.figureNo=99;
-    
-    % ignore zero elements in signal
-    signalRMS=mean(savedInputSignal(find(savedInputSignal)).^2)^0.5;
-    signalRMSdb=20*log10(signalRMS/20e-6);
 
     % plot signal (1)
     plotInstructions.displaydt=dt;
     plotInstructions.numPlots=6;
     plotInstructions.subPlotNo=1;
     plotInstructions.title=...
         ['stimulus (Pa).  ' num2str(signalRMSdb, '%4.0f') ' dB SPL'];
-    r=size(savedInputSignal,1);
-    if r==1, savedInputSignal=savedInputSignal'; end
-    UTIL_plotMatrix(savedInputSignal', plotInstructions);
+    UTIL_plotMatrix(savedInputSignal, plotInstructions);
 
     % stapes (2)
     plotInstructions.subPlotNo=2;
     plotInstructions.title= ['stapes displacement (m)'];
     UTIL_plotMatrix(OMEoutput, plotInstructions);

     plotInstructions.subPlotNo=3;
     plotInstructions.yValues= savedBFlist;
     [r c]=size(DRNLoutput);
     if r>1
         plotInstructions.title= ['BM displacement'];
-    UTIL_plotMatrix(abs(DRNLoutput), plotInstructions);
+        UTIL_plotMatrix(abs(DRNLoutput), plotInstructions);
     else
         plotInstructions.title= ['BM displacement. BF=' ...
             num2str(savedBFlist) ' Hz'];
-    UTIL_plotMatrix(DRNLoutput, plotInstructions);
+        UTIL_plotMatrix(DRNLoutput, plotInstructions);
     end
 
     switch saveAN_spikesOrProbability
         case 'spikes'
             % AN (4)
-            plotInstructions.displaydt=ANdt;
+            plotInstructions.displaydt=dtSpikes;
             plotInstructions.title='AN';
             plotInstructions.subPlotNo=4;
             plotInstructions.yLabel='BF';
             plotInstructions.yValues= savedBFlist;
             plotInstructions.rasterDotSize=1;

                 plotInstructions.rasterDotSize=3;
             end
             UTIL_plotMatrix(ANoutput, plotInstructions);
 
             % CN (5)
-            plotInstructions.displaydt=ANdt;
+            plotInstructions.displaydt=dtSpikes;
             plotInstructions.subPlotNo=5;
             plotInstructions.title='CN spikes';
             if sum(sum(CNoutput))<100
                 plotInstructions.rasterDotSize=3;
             end
             UTIL_plotMatrix(CNoutput, plotInstructions);
 
             % IC (6)
-            plotInstructions.displaydt=ANdt;
+            plotInstructions.displaydt=dtSpikes;
             plotInstructions.subPlotNo=6;
             plotInstructions.title='Brainstem 2nd level';
             if size(ICoutput,1)>1
                 if sum(sum(ICoutput))<100
                     plotInstructions.rasterDotSize=3;

                 plotInstructions.zValuesRange= [-.1 0];
                 UTIL_plotMatrix(ICmembraneOutput, plotInstructions);
             end
 
         otherwise % AN rate based on probability of firing
-            PSTHbinWidth=0.001;
+            PSTHbinWidth=0.0005;
             PSTH= UTIL_PSTHmakerb(ANprobRateOutput, dt, PSTHbinWidth);
-%             PSTH = makeANsmooth(PSTH, 1/PSTHbinWidth);
+            plotInstructions=[];
             plotInstructions.displaydt=PSTHbinWidth;
+            plotInstructions.plotDivider=0;
             plotInstructions.numPlots=2;
             plotInstructions.subPlotNo=2;
             plotInstructions.yLabel='BF';
-            plotInstructions.xLabel='time';
+            plotInstructions.yValues=savedBFlist;
+            plotInstructions.xLabel='time (s)';
             plotInstructions.zValuesRange= [0 300];
+
             if nANfiberTypes>1,
-                plotInstructions.yLabel='LSR    HSR';
+                plotInstructions.yLabel='LSR                   HSR';
                 plotInstructions.plotDivider=1;
             end
             plotInstructions.title='AN - spike rate';
             UTIL_plotMatrix(PSTH, plotInstructions);
             shading interp
-            colorbar('southOutside')
+            if showMapOptions.colorbar, colorbar('southOutside'), end
     end
     set(gcf,'name','MAP output')
 end
 
-if showMapOptions.surfProbability &&...
+
+%% surface plot of AN response
+%   probability
+if showMapOptions.surfAN &&...
         strcmp(saveAN_spikesOrProbability,'probability') && ...
         length(savedBFlist)>2
-    %% surface plot of probability
-        % select only HSR fibers
-        figure(97), clf
-        HSRprobOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:);
-        PSTHbinWidth=0.001;
-        PSTH=UTIL_PSTHmakerb(HSRprobOutput, ANdt, PSTHbinWidth);
-        [nY nX]=size(PSTH);
-        time=PSTHbinWidth*(1:nX);
-        surf(time, savedBFlist, PSTH)
-        shading interp
-        set(gca, 'yScale','log')
-        xlim([0 max(time)])
-        ylim([0 max(savedBFlist)])
-        zlim([0 1000])
-        xlabel('time (s)')
-        ylabel('best frequency (Hz)')
-        zlabel('spike rate')
-        view([-20 60])
-        %     view([0 90])
-        disp(['max max AN: ' num2str(max(max(PSTH)))])
-
-        title (['firing probability of HSR fibers only. Level= ' ...
-            num2str(signalRMSdb,'% 3.0f') ' dB'])
-        colorbar('southOutside')
-end
-
-%% surface plot of AN spikes
-if showMapOptions.surfSpikes ...
-    && strcmp(saveAN_spikesOrProbability, 'spikes')
+    % select only HSR fibers
     figure(97), clf
+    PSTHbinWidth=showMapOptions.PSTHbinwidth;
+    PSTH= UTIL_PSTHmakerb(...
+        ANprobRateOutput(end-length(savedBFlist)+1:end,:), ...
+        dt, PSTHbinWidth);
+    [nY nX]=size(PSTH);
+    time=PSTHbinWidth*(1:nX);
+    surf(time, savedBFlist, PSTH)
+    caxis([0 500])
+    shading interp
+    set(gca, 'yScale','log')
+    xlim([0 max(time)])
+    ylim([0 max(savedBFlist)])
+    zlim([0 500])
+    myFontSize=10;
+    xlabel('time (s)', 'fontsize',myFontSize)
+    ylabel('BF (Hz)', 'fontsize',myFontSize)
+    zlabel('spike rate')
+    view(showMapOptions.view)
+    set(gca,'ytick',[500 1000 2000 8000],'fontSize',myFontSize)
+    title (['AN response. Level= ' ...
+        num2str(signalRMSdb,'% 3.0f') ' dB'...
+        '  binwidth= ' num2str(1000*PSTHbinWidth) ' s'])
+    if showMapOptions.colorbar, colorbar('southOutside'), end
+end
+
+%% surfAN ('spikes')
+if showMapOptions.surfAN ...
+        && strcmp(saveAN_spikesOrProbability, 'spikes')...
+        && length(savedBFlist)>2
+    figure(97), clf
+    % combine fibers across channels
+    nFibersPerChannel=AN_IHCsynapseParams.numFibers;
+    [r nEpochs]=size(ANoutput);
+    nChannels=round(r/nFibersPerChannel);
+    x=reshape(ANoutput,nChannels,nFibersPerChannel,nEpochs);
+    x=squeeze(sum(x,2));
     % select only HSR fibers at the bottom of the matrix
-    ANoutput= ANoutput(end-length(savedBFlist)+1:end,:);
-    PSTHbinWidth=0.005; % 1 ms bins
-    PSTH=UTIL_makePSTH(ANoutput, ANdt, PSTHbinWidth);
+    HSRoutput= x(end-length(savedBFlist)+1:end,:);
+    PSTH=HSRoutput;
+    PSTHbinWidth=showMapOptions.PSTHbinwidth;
+    PSTH=UTIL_makePSTH(HSRoutput, dtSpikes, PSTHbinWidth);
     [nY nX]=size(PSTH);
     time=PSTHbinWidth*(1:nX);
     surf(time, savedBFlist, PSTH)
     shading interp
     set(gca, 'yScale','log')

     ylim([0 max(savedBFlist)])
     %     zlim([0 1000])
     xlabel('time (s)')
     ylabel('best frequency (Hz)')
     zlabel('spike rate')
-    view([-20 60])
-    %     view([0 90])
+    view(showMapOptions.view)
     title ([showMapOptions.fileName ':   ' num2str(signalRMSdb,'% 3.0f') ' dB'])
     set(97,'name', 'spikes surface plot')
 end
 
+%% IC raster plot
+if showMapOptions.ICrasterPlot &&...
+        strcmp(saveAN_spikesOrProbability,'spikes') && ...
+        length(savedBFlist)>2
+    figure(91), clf
+    plotInstructions=[];
+    plotInstructions.numPlots=2;
+    plotInstructions.subPlotNo=2;
+    plotInstructions.title=...
+        ['IC raster plot'];
+    plotInstructions.figureNo=91;
+    plotInstructions.displaydt=dtSpikes;
+    plotInstructions.title='Brainstem 2nd level';
+    plotInstructions.yLabel='BF';
+    plotInstructions.yValues= savedBFlist;
+
+    if size(ICoutput,1)>1
+        if sum(sum(ICoutput))<100
+            plotInstructions.rasterDotSize=3;
+        end
+        UTIL_plotMatrix(ICoutput, plotInstructions);
+    end
+
+    % plot signal (1)
+    plotInstructions.displaydt=dt;
+    plotInstructions.subPlotNo=1;
+    plotInstructions.title=...
+        ['stimulus (Pa).  ' num2str(signalRMSdb, '%4.0f') ' dB SPL'];
+    UTIL_plotMatrix(savedInputSignal, plotInstructions);
+
+end
 
 %% figure(98) plot efferent control values as dB
 if showMapOptions.showEfferent
     plotInstructions=[];
     plotInstructions.figureNo=98;

     plotInstructions.numPlots=4;
     plotInstructions.subPlotNo=1;
     plotInstructions.zValuesRange= [-1 1];
     plotInstructions.title= ['RMS level='...
         num2str(signalRMSdb, '%4.0f') ' dB SPL'];
-    UTIL_plotMatrix(savedInputSignal', plotInstructions);
-    
-    
+    UTIL_plotMatrix(savedInputSignal, plotInstructions);
+
+
     plotInstructions.subPlotNo=2;
     plotInstructions.zValuesRange=[ -25 0];
     plotInstructions.title= ['AR stapes attenuation (dB); tau='...
         num2str(OMEParams.ARtau, '%4.3f') ' s'];
     UTIL_plotMatrix(20*log10(ARattenuation), plotInstructions);

     plotInstructions.subPlotNo=2;
     plotInstructions.yValues= savedBFlist;
     plotInstructions.yLabel= 'dB';
     if strcmp(saveAN_spikesOrProbability,'spikes')
         rate2atten=DRNLParams.rateToAttenuationFactor;
-    plotInstructions.title= ['MOC atten; tau=' ...
-        num2str(DRNLParams.MOCtau) '; factor='...
-        num2str(rate2atten, '%6.4f')];
+        plotInstructions.title= ['MOC atten; tau=' ...
+            num2str(DRNLParams.MOCtau) '; factor='...
+            num2str(rate2atten, '%6.4f')];
     else
         rate2atten=DRNLParams.rateToAttenuationFactorProb;
-    plotInstructions.title= ['MOC atten; tauProb=' ...
-        num2str(DRNLParams.MOCtauProb) '; factor='...
-        num2str(rate2atten, '%6.4f')];
+        plotInstructions.title= ['MOC atten; tauProb=' ...
+            num2str(DRNLParams.MOCtauProb) '; factor='...
+            num2str(rate2atten, '%6.4f')];
     end
     plotInstructions.zValuesRange=[0 -DRNLParams.minMOCattenuationdB+5];
     UTIL_plotMatrix(-20*log10(MOCattenuation), plotInstructions);
     hold on
     [r c]=size(MOCattenuation);
-    if r>2
-    colorbar('southOutside')
+    if r>2 && showMapOptions.colorbar
+        colorbar('southOutside')
     end
     set(plotInstructions.figureNo, 'name', 'AR/ MOC')
-    
+
     binwidth=0.1;
     [PSTH ]=UTIL_PSTHmaker(20*log10(MOCattenuation), dt, binwidth);
     PSTH=PSTH*length(PSTH)/length(MOCattenuation);
     t=binwidth:binwidth:binwidth*length(PSTH);
     fprintf('\n\n')
-%     UTIL_printTabTable([t' PSTH'])
-%     fprintf('\n\n')
-    
-end
-
-%% ACF plot if required
+    %     UTIL_printTabTable([t' PSTH'])
+    %     fprintf('\n\n')
+
+end
+
+%% ACF plot
 if showMapOptions.showACF
     tic
-    method.dt=dt;
-    method.segmentNo=1;
-    method.nonlinCF=savedBFlist;
-
-    minPitch=	80; maxPitch=	4000; numPitches=100;    % specify lags
-    pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches);
-    pitches=fliplr(pitches);
-    filteredSACFParams.lags=1./pitches;     % autocorrelation lags vector
-    filteredSACFParams.acfTau=	.003;       % time constant of running ACF
-    filteredSACFParams.lambda=	0.12;       % slower filter to smooth ACF
-    filteredSACFParams.lambda=	0.01;       % slower filter to smooth ACF
-
-    filteredSACFParams.plotACFs=0;          % special plot (see code)
-    filteredSACFParams.plotFilteredSACF=0;  % 0 plots unfiltered ACFs
-    filteredSACFParams.plotMoviePauses=.3;          % special plot (see code)
-
-    filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags
-    filteredSACFParams.lagsProcedure=  'useAllLags';
-    % filteredSACFParams.lagsProcedure=  'useBernsteinLagWeights';
-    % filteredSACFParams.lagsProcedure=  'omitShortLags';
-    filteredSACFParams.criterionForOmittingLags=3;
-    filteredSACFParams.plotACFsInterval=200;
-
     if filteredSACFParams.plotACFs
         % plot original waveform on ACF plot
         figure(13), clf
         subplot(4,1,1)
         t=dt*(1:length(savedInputSignal));
         plot(t,savedInputSignal)
         xlim([0 t(end)])
         title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']);
     end
 
+    % compute ACF
+    switch saveAN_spikesOrProbability
+        case 'probability'
+            inputToACF=ANprobRateOutput(end-length(savedBFlist)+1:end,:);
+        otherwise
+            inputToACF=ANoutput;
+    end
+
+    disp ('computing ACF...')
+    [P, BFlist, sacf, boundaryValue] = ...
+        filteredSACF(inputToACF, dt, savedBFlist, filteredSACFParams);
+    disp(' ACF done.')
+
     % plot original waveform on summary/smoothed ACF plot
     figure(96), clf
     subplot(2,1,1)
     t=dt*(1:length(savedInputSignal));
     plot(t,savedInputSignal)
     xlim([0 t(end)])
     title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']);
 
-
-    % compute ACF
-    switch saveAN_spikesOrProbability
-        case 'probability'
-            inputToACF=ANprobRateOutput.^0.5;
-        otherwise
-            inputToACF=ANoutput;
-    end
-
-    disp ('computing ACF...')
-    [P, BFlist, sacf, boundaryValue] = ...
-        filteredSACF(inputToACF, method, filteredSACFParams);
-    disp(' ACF done.')
-
-    % SACF
+    % plot SACF
+    figure(96)
     subplot(2,1,2)
     imagesc(P)
+%     surf(filteredSACFParams.lags, t, P)
     ylabel('periodicities (Hz)')
     xlabel('time (s)')
-    title(['running smoothed (root) SACF. ' saveAN_spikesOrProbability ' input'])
-    pt=[1 get(gca,'ytick')]; % force top xtick to show
+    title(['running smoothed SACF. ' saveAN_spikesOrProbability ' input'])
+    pt=[1 get(gca,'ytick')]; % force top ytick to show
     set(gca,'ytick',pt)
+    pitches=1./filteredSACFParams.lags;     % autocorrelation lags vector
     set(gca,'ytickLabel', round(pitches(pt)))
+    [nCH nTimes]=size(P);
+    t=dt:dt:dt*nTimes;
     tt=get(gca,'xtick');
     set(gca,'xtickLabel', round(100*t(tt))/100)
 end
 
 path(restorePath)

 % figure(95), plot(CNtauGk,ICrate)
 % title ('ICrate'), xlabel('CNtauGk'), ylabel('ICrate')
 % end
 
 function ANsmooth = makeANsmooth(ANresponse, sampleRate, winSize, hopSize)
-            if nargin < 3
-                winSize = 25; %default 25 ms window
-            end
-            if nargin < 4
-                hopSize = 10; %default 10 ms jump between windows
-            end
-            
-            winSizeSamples = round(winSize*sampleRate/1000);
-            hopSizeSamples = round(hopSize*sampleRate/1000);
-            
-            % smooth
-            hann = hanning(winSizeSamples);
-            
-            ANsmooth = [];%Cannot pre-allocate a size as it is unknown until the enframing
-            for chan = 1:size(ANresponse,1)
-                f = enframe(ANresponse(chan,:), hann, hopSizeSamples);
-                ANsmooth(chan,:) = mean(f,2)'; %#ok<AGROW> see above comment
-            end
+if nargin < 3
+    winSize = 25; %default 25 ms window
+end
+if nargin < 4
+    hopSize = 10; %default 10 ms jump between windows
+end
+
+winSizeSamples = round(winSize*sampleRate/1000);
+hopSizeSamples = round(hopSize*sampleRate/1000);
+
+% smooth
+hann = hanning(winSizeSamples);
+
+ANsmooth = [];%Cannot pre-allocate a size as it is unknown until the enframing
+for chan = 1:size(ANresponse,1)
+    f = enframe(ANresponse(chan,:), hann, hopSizeSamples);
+    ANsmooth(chan,:) = mean(f,2)'; %#ok<AGROW> see above comment
+end
 %         end% ------ OF makeANsmooth