rmeddis@26: function UTIL_showMAP (showMapOptions, paramChanges) rmeddis@24: % UTIL_showMAP produces summaries of the output from MAP's mmost recent run rmeddis@24: % All MAP outputs are stored in global variables and UTIL_showMAP rmeddis@24: % simply assumes that they are in place. rmeddis@24: % rmeddis@25: % showMapOptions rmeddis@25: % showMapOptions.printModelParameters=1; % print model parameters rmeddis@25: % showMapOptions.showModelOutput=1; % plot all stages output rmeddis@25: % showMapOptions.printFiringRates=1; % mean activity at all stages rmeddis@25: % showMapOptions.showACF=1; % SACF (probabilities only) rmeddis@25: % showMapOptions.showEfferent=1; % plot of efferent activity rmeddis@25: % showMapOptions.surfProbability=0; % HSR (probability) surf plot rmeddis@25: % showMapOptions.fileName=[]; % parameter filename for plot title rmeddis@23: rmeddis@23: dbstop if warning rmeddis@23: rmeddis@23: global dt ANdt saveAN_spikesOrProbability savedBFlist saveMAPparamsName... rmeddis@23: savedInputSignal TMoutput OMEoutput ARattenuation ... rmeddis@23: DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV... rmeddis@23: IHCoutput ANprobRateOutput ANoutput savePavailable tauCas ... rmeddis@23: CNoutput ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation rmeddis@23: global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams rmeddis@23: global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams rmeddis@23: rmeddis@23: rmeddis@23: restorePath=path; rmeddis@23: addpath ( ['..' filesep 'utilities'], ['..' filesep 'parameterStore']) rmeddis@23: rmeddis@23: if nargin<1 rmeddis@25: showMapOptions=[]; rmeddis@23: end rmeddis@25: % defaults (plot staged outputs and print rates only) if options omitted rmeddis@25: if ~isfield(showMapOptions,'printModelParameters') rmeddis@25: showMapOptions.printModelParameters=0; end rmeddis@25: if ~isfield(showMapOptions,'showModelOutput'),showMapOptions.showModelOutput=1;end rmeddis@25: if ~isfield(showMapOptions,'printFiringRates'),showMapOptions.printFiringRates=1;end rmeddis@25: if ~isfield(showMapOptions,'showACF'),showMapOptions.showACF=0;end rmeddis@25: if ~isfield(showMapOptions,'showEfferent'),showMapOptions.showEfferent=0;end rmeddis@25: if ~isfield(showMapOptions,'surfProbability'),showMapOptions.surfProbability=0;end rmeddis@25: if ~isfield(showMapOptions,'fileName'),showMapOptions.fileName=[];end rmeddis@25: if ~isfield(showMapOptions,'surfSpikes'),showMapOptions.surfSpikes=0;end rmeddis@23: rmeddis@26: %% send all model parameters to command window rmeddis@25: if showMapOptions.printModelParameters rmeddis@23: % Read parameters from MAPparams<***> file in 'parameterStore' folder rmeddis@23: % and print out all parameters rmeddis@26: if nargin>1 rmeddis@26: cmd=['MAPparams' saveMAPparamsName '(-1, 1/dt, 1, paramChanges);']; rmeddis@26: eval(cmd); rmeddis@26: else rmeddis@26: cmd=['MAPparams' saveMAPparamsName '(-1, 1/dt, 1);']; rmeddis@26: eval(cmd); rmeddis@26: disp(' no parameter changes found') rmeddis@26: end rmeddis@23: end rmeddis@23: rmeddis@26: %% summarise firing rates in command window rmeddis@25: if showMapOptions.printFiringRates rmeddis@23: %% print summary firing rates rmeddis@23: fprintf('\n\n') rmeddis@23: disp('summary') rmeddis@23: disp(['AR: ' num2str(min(ARattenuation))]) rmeddis@23: disp(['MOC: ' num2str(min(min(MOCattenuation)))]) rmeddis@23: nANfiberTypes=length(tauCas); rmeddis@23: if strcmp(saveAN_spikesOrProbability, 'spikes') rmeddis@23: nANfibers=size(ANoutput,1); rmeddis@23: nHSRfibers=nANfibers/nANfiberTypes; rmeddis@23: duration=size(TMoutput,2)*dt; rmeddis@23: disp(['AN: ' num2str(sum(sum(ANoutput(end-nHSRfibers+1:end,:)))/... rmeddis@23: (nHSRfibers*duration))]) rmeddis@23: rmeddis@23: nCNneurons=size(CNoutput,1); rmeddis@23: nHSRCNneuronss=nCNneurons/nANfiberTypes; rmeddis@23: disp(['CN: ' num2str(sum(sum(CNoutput(end-nHSRCNneuronss+1:end,:)))... rmeddis@23: /(nHSRCNneuronss*duration))]) rmeddis@23: disp(['IC: ' num2str(sum(sum(ICoutput))/duration)]) rmeddis@23: % disp(['IC by type: ' num2str(mean(ICfiberTypeRates,2)')]) rmeddis@23: else rmeddis@27: HSRprobOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:); rmeddis@27: disp(['AN(HSR): ' num2str(mean(mean(HSRprobOutput)))]) rmeddis@27: PSTH= UTIL_PSTHmakerb(HSRprobOutput, dt, 0.001); rmeddis@27: disp(['max max AN: ' num2str(max(max(PSTH)))]) rmeddis@23: end rmeddis@23: end rmeddis@23: rmeddis@23: rmeddis@26: %% figure (99): display output from all model stages rmeddis@25: if showMapOptions.showModelOutput rmeddis@23: plotInstructions.figureNo=99; rmeddis@23: signalRMS=mean(savedInputSignal.^2)^0.5; rmeddis@23: signalRMSdb=20*log10(signalRMS/20e-6); rmeddis@23: rmeddis@23: % plot signal (1) rmeddis@23: plotInstructions.displaydt=dt; rmeddis@23: plotInstructions.numPlots=6; rmeddis@23: plotInstructions.subPlotNo=1; rmeddis@23: plotInstructions.title=... rmeddis@23: ['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']; rmeddis@23: r=size(savedInputSignal,1); rmeddis@23: if r==1, savedInputSignal=savedInputSignal'; end rmeddis@23: UTIL_plotMatrix(savedInputSignal', plotInstructions); rmeddis@23: rmeddis@23: % stapes (2) rmeddis@23: plotInstructions.subPlotNo=2; rmeddis@23: plotInstructions.title= ['stapes displacement']; rmeddis@23: UTIL_plotMatrix(OMEoutput, plotInstructions); rmeddis@23: rmeddis@23: % DRNL (3) rmeddis@23: plotInstructions.subPlotNo=3; rmeddis@23: plotInstructions.title= ['BM displacement']; rmeddis@23: plotInstructions.yValues= savedBFlist; rmeddis@23: UTIL_plotMatrix(DRNLoutput, plotInstructions); rmeddis@23: rmeddis@23: switch saveAN_spikesOrProbability rmeddis@23: case 'spikes' rmeddis@23: % AN (4) rmeddis@23: plotInstructions.displaydt=ANdt; rmeddis@23: plotInstructions.title='AN'; rmeddis@23: plotInstructions.subPlotNo=4; rmeddis@23: plotInstructions.yLabel='BF'; rmeddis@23: plotInstructions.yValues= savedBFlist; rmeddis@23: plotInstructions.rasterDotSize=1; rmeddis@23: if length(tauCas)==2 rmeddis@23: plotInstructions.plotDivider=1; rmeddis@23: else rmeddis@23: plotInstructions.plotDivider=0; rmeddis@23: end rmeddis@23: if sum(sum(ANoutput))<100 rmeddis@23: plotInstructions.rasterDotSize=3; rmeddis@23: end rmeddis@23: UTIL_plotMatrix(ANoutput, plotInstructions); rmeddis@23: rmeddis@23: % CN (5) rmeddis@23: plotInstructions.displaydt=ANdt; rmeddis@23: plotInstructions.subPlotNo=5; rmeddis@23: plotInstructions.title='CN spikes'; rmeddis@23: if sum(sum(CNoutput))<100 rmeddis@23: plotInstructions.rasterDotSize=3; rmeddis@23: end rmeddis@23: UTIL_plotMatrix(CNoutput, plotInstructions); rmeddis@23: rmeddis@23: % IC (6) rmeddis@23: plotInstructions.displaydt=ANdt; rmeddis@23: plotInstructions.subPlotNo=6; rmeddis@23: plotInstructions.title='IC'; rmeddis@23: if size(ICoutput,1)>3 rmeddis@23: if sum(sum(ICoutput))<100 rmeddis@23: plotInstructions.rasterDotSize=3; rmeddis@23: end rmeddis@23: UTIL_plotMatrix(ICoutput, plotInstructions); rmeddis@23: else rmeddis@23: plotInstructions.title='IC (HSR) membrane potential'; rmeddis@23: plotInstructions.displaydt=dt; rmeddis@23: plotInstructions.yLabel='V'; rmeddis@23: plotInstructions.zValuesRange= [-.1 0]; rmeddis@23: UTIL_plotMatrix(ICmembraneOutput, plotInstructions); rmeddis@23: end rmeddis@23: rmeddis@26: otherwise % AN rate based on probability of firing rmeddis@26: PSTHbinWidth=0.001; rmeddis@26: PSTH= UTIL_PSTHmakerb(ANprobRateOutput, dt, PSTHbinWidth); rmeddis@26: plotInstructions.displaydt=PSTHbinWidth; rmeddis@23: plotInstructions.numPlots=2; rmeddis@23: plotInstructions.subPlotNo=2; rmeddis@23: plotInstructions.yLabel='BF'; rmeddis@23: if nANfiberTypes>1, rmeddis@23: plotInstructions.yLabel='LSR HSR'; rmeddis@23: plotInstructions.plotDivider=1; rmeddis@23: end rmeddis@26: plotInstructions.title='AN - spike rate'; rmeddis@26: UTIL_plotMatrix(PSTH, plotInstructions); rmeddis@23: end rmeddis@23: end rmeddis@23: rmeddis@25: if showMapOptions.surfProbability rmeddis@23: %% surface plot of probability rmeddis@26: if strcmp(saveAN_spikesOrProbability,'probability') && ... rmeddis@26: length(savedBFlist)>2 rmeddis@23: figure(97), clf rmeddis@23: % select only HSR fibers at the bottom of the matrix rmeddis@27: HSRprobOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:); rmeddis@23: PSTHbinWidth=0.001; rmeddis@27: PSTH=UTIL_PSTHmakerb(HSRprobOutput, ANdt, PSTHbinWidth); rmeddis@23: [nY nX]=size(PSTH); rmeddis@23: time=PSTHbinWidth*(1:nX); rmeddis@23: surf(time, savedBFlist, PSTH) rmeddis@23: shading interp rmeddis@23: set(gca, 'yScale','log') rmeddis@23: xlim([0 max(time)]) rmeddis@23: ylim([0 max(savedBFlist)]) rmeddis@27: zlim([0 2000]) rmeddis@23: xlabel('time (s)') rmeddis@23: ylabel('best frequency (Hz)') rmeddis@23: zlabel('spike rate') rmeddis@23: view([-20 60]) rmeddis@23: % view([0 90]) rmeddis@27: disp(['max max AN: ' num2str(max(max(PSTH)))]) rmeddis@27: rmeddis@25: title ([showMapOptions.fileName ': ' num2str(signalRMSdb,'% 3.0f') ' dB']) rmeddis@26: end rmeddis@23: end rmeddis@23: rmeddis@25: if showMapOptions.surfSpikes rmeddis@26: %% surface plot of AN spikes rmeddis@25: figure(97), clf rmeddis@25: % select only HSR fibers at the bottom of the matrix rmeddis@25: ANoutput= ANoutput(end-length(savedBFlist)+1:end,:); rmeddis@26: PSTHbinWidth=0.005; % 1 ms bins rmeddis@26: PSTH=UTIL_makePSTH(ANoutput, ANdt, PSTHbinWidth); rmeddis@25: [nY nX]=size(PSTH); rmeddis@25: time=PSTHbinWidth*(1:nX); rmeddis@25: surf(time, savedBFlist, PSTH) rmeddis@25: shading interp rmeddis@25: set(gca, 'yScale','log') rmeddis@25: xlim([0 max(time)]) rmeddis@25: ylim([0 max(savedBFlist)]) rmeddis@25: % zlim([0 1000]) rmeddis@25: xlabel('time (s)') rmeddis@25: ylabel('best frequency (Hz)') rmeddis@25: zlabel('spike rate') rmeddis@25: view([-20 60]) rmeddis@25: % view([0 90]) rmeddis@25: title ([showMapOptions.fileName ': ' num2str(signalRMSdb,'% 3.0f') ' dB']) rmeddis@25: end rmeddis@25: rmeddis@23: rmeddis@26: %% figure(98) plot efferent control values as dB rmeddis@25: if showMapOptions.showEfferent rmeddis@23: plotInstructions=[]; rmeddis@23: plotInstructions.figureNo=98; rmeddis@23: figure(98), clf rmeddis@23: plotInstructions.displaydt=dt; rmeddis@23: plotInstructions.numPlots=2; rmeddis@23: plotInstructions.subPlotNo=1; rmeddis@23: plotInstructions.zValuesRange=[ -25 0]; rmeddis@23: plotInstructions.title= ['AR strength. Signal level= ' ... rmeddis@23: num2str(signalRMSdb,'%4.0f') ' dB SPL']; rmeddis@23: UTIL_plotMatrix(20*log10(ARattenuation), plotInstructions); rmeddis@23: rmeddis@23: plotInstructions.subPlotNo=2; rmeddis@23: plotInstructions.yValues= savedBFlist; rmeddis@23: plotInstructions.yLabel= 'BF'; rmeddis@23: plotInstructions.title= ['MOC strength']; rmeddis@23: plotInstructions.zValuesRange=[ -25 0]; rmeddis@23: subplot(2,1,2) rmeddis@23: % imagesc(MOCattenuation) rmeddis@23: UTIL_plotMatrix(20*log10(MOCattenuation), plotInstructions); rmeddis@23: colorbar rmeddis@23: end rmeddis@23: rmeddis@23: %% ACF plot if required rmeddis@25: if showMapOptions.showACF rmeddis@23: tic rmeddis@23: method.dt=dt; rmeddis@23: method.segmentNo=1; rmeddis@23: method.nonlinCF=savedBFlist; rmeddis@23: rmeddis@23: minPitch= 80; maxPitch= 4000; numPitches=100; % specify lags rmeddis@23: pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches); rmeddis@23: pitches=fliplr(pitches); rmeddis@23: filteredSACFParams.lags=1./pitches; % autocorrelation lags vector rmeddis@23: filteredSACFParams.acfTau= .003; % time constant of running ACF rmeddis@23: filteredSACFParams.lambda= 0.12; % slower filter to smooth ACF rmeddis@23: filteredSACFParams.lambda= 0.01; % slower filter to smooth ACF rmeddis@23: rmeddis@23: filteredSACFParams.plotACFs=0; % special plot (see code) rmeddis@23: filteredSACFParams.plotFilteredSACF=0; % 0 plots unfiltered ACFs rmeddis@23: filteredSACFParams.plotMoviePauses=.3; % special plot (see code) rmeddis@23: rmeddis@23: filteredSACFParams.usePressnitzer=0; % attenuates ACF at long lags rmeddis@23: filteredSACFParams.lagsProcedure= 'useAllLags'; rmeddis@23: % filteredSACFParams.lagsProcedure= 'useBernsteinLagWeights'; rmeddis@23: % filteredSACFParams.lagsProcedure= 'omitShortLags'; rmeddis@23: filteredSACFParams.criterionForOmittingLags=3; rmeddis@23: filteredSACFParams.plotACFsInterval=200; rmeddis@23: rmeddis@23: if filteredSACFParams.plotACFs rmeddis@23: % plot original waveform on ACF plot rmeddis@23: figure(13), clf rmeddis@23: subplot(4,1,1) rmeddis@23: t=dt*(1:length(savedInputSignal)); rmeddis@23: plot(t,savedInputSignal) rmeddis@23: xlim([0 t(end)]) rmeddis@23: title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); rmeddis@23: end rmeddis@23: rmeddis@23: % plot original waveform on summary/smoothed ACF plot rmeddis@23: figure(96), clf rmeddis@23: subplot(2,1,1) rmeddis@23: t=dt*(1:length(savedInputSignal)); rmeddis@23: plot(t,savedInputSignal) rmeddis@23: xlim([0 t(end)]) rmeddis@23: title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); rmeddis@23: rmeddis@23: rmeddis@23: % compute ACF rmeddis@23: switch saveAN_spikesOrProbability rmeddis@23: case 'probability' rmeddis@23: inputToACF=ANprobRateOutput.^0.5; rmeddis@23: otherwise rmeddis@23: inputToACF=ANoutput; rmeddis@23: end rmeddis@23: rmeddis@23: disp ('computing ACF...') rmeddis@23: [P, BFlist, sacf, boundaryValue] = ... rmeddis@23: filteredSACF(inputToACF, method, filteredSACFParams); rmeddis@23: disp(' ACF done.') rmeddis@23: rmeddis@23: % SACF rmeddis@23: subplot(2,1,2) rmeddis@23: imagesc(P) rmeddis@23: ylabel('periodicities (Hz)') rmeddis@23: xlabel('time (s)') rmeddis@23: title(['running smoothed (root) SACF. ' saveAN_spikesOrProbability ' input']) rmeddis@23: pt=[1 get(gca,'ytick')]; % force top xtick to show rmeddis@23: set(gca,'ytick',pt) rmeddis@23: set(gca,'ytickLabel', round(pitches(pt))) rmeddis@23: tt=get(gca,'xtick'); rmeddis@23: set(gca,'xtickLabel', round(100*t(tt))/100) rmeddis@23: end rmeddis@23: rmeddis@23: path(restorePath)