function neuroplot( varargin )

%NEUROPLOT  Plot instantaneous or window-averaged Vm and spike firing rate patterns over time. 
%           Create .avi movies and/or .wav sounds recording of spike patterns as they evolve.
%           Concatenates multiple wildcard-matching *.txt files into a single .mat, saved for later upload.
%
% rev 30June2002 Phil Goodman    goodman@unr.edu    brain.unr.edu 
% 30June2002: changed default spikepeak to 30 mV; fixed raw2d colorbar legend scaling
%
% Use the following 'parameter', 'value' pairs in any order...
%
% fname         'filename.txt' or 'filename.mat' with matrix of cell membrane voltages
%                   File format: Rectangular data, separated by space or tab (see also celldata, below)
%                                First column (row, if 'celldata', 'row') may be serial integers labeling tick
%                                (will be ignored). 
%                   'filename.mat' is saved by neuroplot() for later, faster, uploading; cells are always in
%                                columns, and any preexisting tick column label is stripped.
%
% OPTIONAL ARGUMENTS
% celldata      'columns' or 'rows', indicates that columns (rows) correspond to cells [DEFAULT 'columns']
%                  Will ignore first column (row) if it contains sequential integers, assuming it's a tick label.
% savemat       'yes' or 'no' to save for later quick load as  fname.mat [only if fname is of type .txt] 
%                  (Multiple file concatenation: if xxx*.txt used, * will replaced by _ )  [DEFAULT 'yes']
% plotmode      'raw2d' displays raw data in pseudocolor matrix
%               'raw3d' displays raw data in 3d landscape
%               'rawdot' displays spikes at dots in 2d raster matrix
%               'rawsubplots' displays each cell Vm tracing in a separate subplot
%               'avg2d' displays window-averaged firing rates in pseudocolor matrix
%               'avg3d' displays window-averaged firing in 3d landscape
%               'density' displays pseudocolored contour of spike density (uses 'tick:cell_compression')
% dotsize       [integer] size of dot in rawdot plot [ DEFAULT 5 ]
% stuntspike    'no' or 'yes' to shrink spike height (raw3d only) to optimize view of subthresh activity [DEFAULT 'no' ]
% viewazel      [integer, integer] for 3d objects azimuth and elevation of view [DEFAULT [70,45] ]
% cmap          'jet', 'hot' 'hsv', 'cool', etc for colormap of color figures [DEFAULT 'jet', except 'hot' for avg2d ]
% tickcompress  (integer) compression of ticks, only used for density plot [DEFAULT 50]
% cellcompress  (integer) compression of cells, only used for density plot [DEFAULT 50]
% window        (integer) window, in ticks, to average (avg* modes only) [DEFAULT 0.1sec * fsv]
% cellrange     'from:to' or 'from:by:to' range of cells to sample (skipping "by") [DEFAULT '1:end']
%                 or 'idfname', filename containing vector or string of cell id number to use.
% tickrange     'from:to' or 'from:by:to' range of ticks to sample (skipping "by") [DEFAULT '1:end']
% refrange      'from:to', reference range averaged, subtracted from comprange to sort 'avg*' plot_types
% comprange     'from:to', reference range averaged, for 'avg*' plots (uses  averages alone if no 'refrange' specified)
%                 (Note: If refrange and comprange not defined, no sorting will be done)
% fsv           (integer) sampling freq, in ticks/second [DEFAULT 10,000 ]
% thresh        (float) threshold for spike, mV [DEFAULT -40]
% spikepeak     (float) peak of action potential, in millivolts [DEFAULT 40.0 ]
% play          'none', 'sound', 'movie', 'both', [DEFAULT 'none']
%                   Notes: 'sound' available only for raw2d, raw3d, and rawdot plotmodes
%                          'movie' available only for raw2d, raw3d, rawdot, avg2d, and avg3d plotmodes
% soundfreq     (integer) for sound, frequency of audio play rate sampling [DEFAULT 3000]
% soundsubsample 'from:by:to' for sound, subsample the plotted cells in this range [DEFAULT '1:end']
%                'random' will randomly select, at each tick, one of the cells for activity for sound
% nframes       (integer) for movies, numbers of frames into which to split figure [DEFAULT 64 ]
% framerate     (integer) for movies, playback at this rate per second [DEFAULT 30 ]
% basecolor     'on' vs 'off' for movies (advance color of plot area not yet filled) [DEFAULT 'on', except 'raw3d']
% xtick2dfactor (integer) by this factor, splits tick into narrower intervals than usual Matlab plot [DEFAULT: 2 ]
% figpos        [ left bottom width height ] (integers) fig position, 0-1 normalized to top left screen [DEFAULT [.1 .5  .8  .35]]
% clearfig      'yes' or 'no'  to clear existing figures at start of this run
% deletenan     'no' or 'yes', to delete cells if any NaNs found within the tickrange [DEFAULT 'no']
%
% EXAMPLES:     neuroplot( 'fname', 'report_01.txt')
%                  Performs raw2d plot on cell data in columns of 'report_01.txt'
%
%               neuroplot( 'fname', 'report_01.txt', 'plotmode', 'raw3d', 'celldata', 'rows', 'play', 'movie' )
%                  Performs raw3d plot on cell data in rows of 'report_01.txt', and creates a movie thereof.

fprintf(1,'\n--> Starting NEUROPLOT (last rev. 30July01; see brain.cs.unr.edu for latest version) ...')
% Assign arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
legal_tokens = { 'fname', 'celldata', 'savemat', 'plotmode',  'window', 'cellrange', 'tickrange','refrange', ...
        'comprange', 'spikepeak', 'fsv', 'play' ,'soundfreq', 'soundsubsample', 'nframes', 'framerate', ...
        'basecolor', 'tickcompress', 'cellcompress', 'viewazel', 'cmap', 'xtick2dfactor','figpos',...
        'deletenan', 'dotsize','clearfig', 'stuntspike', 'thresh' };
[ tokenlist, valuelist ] = parsearg( varargin, legal_tokens );
for t = 1:length(tokenlist), eval([ tokenlist{t} '=valuelist{t};' ]), end 

% DEFAULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~exist( 'celldata' ), celldata = 'columns'; end
if ~exist( 'savemat' ), savemat = 'yes'; end
if ~exist( 'plotmode' ), plotmode = 'raw2d'; end
if ~exist( 'cellrange' ), cellrange = '1:end'; end
if ~exist( 'tickrange' ), tickrange = '1:end'; end
if ~exist( 'refrange' ) & ~exist( 'comprange' ) % neither or both must be defined
    refrange = 'nosort';
    comprange = 'nosort'; 
elseif ~exist( 'refrange' ) & exist( 'comprange' )
       refrange = 'nosort'; % absolute values of comprange, no refrange
end
if ~exist( 'fsv' ), fsv = 1E4; end
if ~exist( 'window' ), window = fsv*0.1; end % 100 ms default window
if ~exist( 'spikepeak' ), spikepeak = 30; end
if ~exist( 'play' ), play = 'none'; end
if ~exist( 'soundfreq' ), soundfreq = 3000; end
if ~exist( 'soundsubsample' ), soundsubsample = '1:end'; end
if ~exist( 'nframes' ), nframes = 64; end
if ~exist( 'basecolor' )
    if strcmp( plotmode, 'raw3d'), basecolor = 'off'; else basecolor = 'on'; end
end
if ~exist( 'framerate' ), framerate = 6; end
if ~exist( 'tickcompress' ), tickcompress = 50; end
if ~exist( 'cellcompress' ), cellcompress = 50; end
if ~exist( 'viewazel' ), viewazel = [70,45]; end
if ~exist( 'cmap' ), cmap = 'jet'; end
if ~exist( 'xtick2dfactor' ), xtick2dfactor = 2; end
if ~exist( 'figpos' ), figpos = [0.1    0.5    0.8    0.35]; end
if ~exist( 'deletenan' ), deletenan = 'no'; end
if ~exist( 'dotsize' ), dotsize = 5; end
if ~exist( 'clearfig' ), clearfig = 'no'; end
if ~exist( 'stuntspike' ), stuntspike = 'no'; end
if ~exist( 'thresh' ), thresh = -40; end

% clear existing figures if requested
if clearfig(1) == 'y'
    cf;
end

% Determine if multiple files to be concatenated %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
starindex = findstr('*', fname );
if strcmp(fname(end-3:end), '.mat') | isempty( starindex )
    flist = {};
else
    flist = get_list_of_fnames(fname, starindex);
end

plot_learning( flist, fname, celldata, savemat, cellrange, tickrange, window,...
    refrange, comprange, spikepeak, fsv, plotmode, play, soundfreq, soundsubsample,nframes, framerate,...
    basecolor, tickcompress, cellcompress, viewazel, cmap, xtick2dfactor, figpos, deletenan, dotsize,...
    stuntspike, thresh );
fprintf(1,'\n--> DONE. Bye!')

% SUBROUTINES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot_learning()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_learning( flist, fname, celldata, savemat, cellrange, tickrange, window,...
    refrange, comprange, spikepeak, fsv, plotmode, play, soundfreq, soundsubsample,nframes, framerate,...
    basecolor, tickcompress, cellcompress, viewazel, cmap, xtick2dfactor, figpos, deletenan, dotsize,...
    stuntspike, thresh );

check_lessthan_minus80 = 1;
if isempty( flist ) % single .txt or .mat file input
    matname = [ fname(1:end-4) '.mat' ];
    extension = fname(end-3:end);
    loadname = fname;
    if strcmp(extension, '.mat')  % already saved as .mat
        fprintf(1,'\n--> Loading EXISTING MATFILE ''%s''...', matname)
        omatrix = load( matname );
        omatrix = struct2cell( omatrix );
        omatrix = omatrix{:};
        if deletenan(1) == 'y' % for previously saved .mat, check/delete nans only if requested
            omatrix = checkfornans( loadname, omatrix, deletenan );       
        end
    else
        fprintf(1,'\n--> Loading NEW TEXT file ''%s''...', loadname)
        omatrix = load( loadname, '-ascii');
        omatrix = test_tick_col(omatrix, celldata);  % strip first column if tick labels
        if strcmp(savemat,'yes')
                save( matname, 'omatrix' ); 
                fprintf(1,'\n--> Saved ''%s'' for future loading, celldata now in columns.', matname )
        end 
        % test for NaNs  (not tested if .mat file is loaded, as was prev reported when .txt loaded)
        omatrix = checkfornans( loadname, omatrix, deletenan );       
    end
    
else % multiple .txt files to concatenate
    % get first file and preallocate memory
    fprintf(1,'\n--> Concatenating MULTIPLE TEXT files using ''%s''...', fname )
    fdata = load( flist{1}, '-ascii');
    fdata = test_tick_col( fdata, celldata );
    fprintf(1,'\n--> Loaded ''%s'', stripped size(%d,%d)...', flist{1}, size(fdata,1), size(fdata,2))
    
    % test for NaNs  (not tested if .mat file is loaded, as was prev reported when .txt loaded)
    fdata = checkfornans( flist{1}, fdata, 'no');       
    nticksperfile = size(fdata, 1);
    ncells = size(fdata, 2);
    nrows = nticksperfile*length(flist);
    %omatrix = repmat( nan, nrows, size(fdata, 2) ); % preallocate memory to speed loading
    omatrix = fdata; % temp dynamic alloc because not sure about each file tick length
    for f = 2:length(flist)
        fdata = load( flist{f}, '-ascii');
        fdata = test_tick_col( fdata, celldata );
        fprintf(1,'\n--> Loaded ''%s'', stripped size(%d,%d)...', flist{f}, size(fdata,1), size(fdata,2))
        if size(fdata, 1) ~= nticksperfile
            fprintf(1,'\n!!! Ticks (%d) in %s is not same as first file (%d ticks).',size(fdata, 1),flist{f},nticksperfile);
        end
        if size(fdata, 2) ~= ncells
            fprintf(1,'\n!!! Cell (%d) in %s is not same as first file (%d cells).',size(fdata, 2),flist{f},ncells);            
        end
        fdata = checkfornans( flist{f}, fdata, 'no' );       
        omatrix = [ omatrix; fdata ];
    end
    matname = [ fname(1:end-5) '~.mat' ]; % replace * with ~
    if strcmp(savemat,'yes')
        save( matname, 'omatrix' ); 
        fprintf(1,'\n--> Saved ''%s'' for future loading.', matname )
    end 
    if deletenan(1) == 'y' % for previously saved .mat, check/delete nans only if requested
        omatrix = checkfornans( loadname, omatrix, deletenan );       
    end
end

pack % compress memory
% select tick range and cells
if isempty( find( cellrange == ':' ) ) % if cellrange is in a file, convert
    cellrangenum = load( cellrange,  '-ascii' );
    cellrangenum = cellrangenum -  nanmin( cellrangenum ) + 1; % reset first index as 1
    omatrix = eval( ['omatrix( ' tickrange ', cellrangenum );'] );
else
    omatrix = eval( ['omatrix( ' tickrange ',' cellrange ');'] );
end

if check_lessthan_minus80 % not for ~.mat files, as these values already replaced with -80 when loaded as .txt
    % limit lower range to -80 and report if so
    lowindices = find( omatrix < -80 );
    if ~isempty( lowindices )
        omatrix( lowindices ) = -80;
        nlowcells = 0;
        for ocol = 1:size(omatrix,2)
            islow = any( omatrix( :, ocol ) < -80 );
            nlowcells = nlowcells + islow;
        end
        fprintf(1,'\n--> NOTE: Found %d cells (w/ %d total ticks) < -80 mV.',...
            nlowcells, length(lowindices) );
    end
end

% for titles, escape titles to prevent underscore
fname = slashunderscore( fname );

% RAW PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(strmatch(plotmode,{'raw2d', 'raw3d'} ) )% do color raw data
    pos = [155         690        1395         420];
    XX = [1:size( omatrix, 1)];
    YY = [1:size( omatrix, 2)];
    xmax = size( omatrix, 1 );
    ymax = size( omatrix, 2 ); 
    zmin = nanmin( omatrix(:)); zmax = nanmax( omatrix(:));
    if strcmp(plotmode,'raw3d')
        % shrink spike height if requested for more optimal viewing subthreshold activity
        if stuntspike(1) == 'y'
            nticks = size(omatrix,1);
            replacement = [5.0000   10.0000   7.0000  2.0000]';
            for mcol = 1:size(omatrix,2)
                sindexes = find( omatrix( :, mcol) == 21);  % lay down shorter spike here
                if ~isempty( sindexes)
                    for si = 1:length( sindexes )
                        starttick = sindexes(si);
                        endtick = starttick + 3;
                        if endtick <= nticks
                            spikeind = [  sindexes(si) : sindexes(si) + 3 ];
                            omatrix( spikeind, mcol) = replacement;
                        else
                            spikeind = [sindexes(si):(sindexes(si) + (3-(endtick-nticks))) ];
                            omatrix( spikeind, mcol) = replacement(1: (4-(endtick-nticks)) );
                        end
                    end
                end
            end
        end
        snfig = figure('units', 'normalized','position', figpos);
        if ~isempty(strmatch(play,{'movie', 'both'} ) )
            fprintf(1,'\n--> Making MOVIE of %d frames to play at %d /sec...', nframes, framerate)
            movfname = [ matname(1:end-4) '.avi' ];
            ticksperframe = ceil( size( omatrix, 1 )/nframes );
            for k = 1:nframes
                endtick = k*ticksperframe;
                if endtick > size( omatrix, 1 )
                    endtick = size( omatrix, 1 );
                end     
                if strcmp(basecolor,'on') % forces a -60 mV blue baseline rather than transparent
                    tempmatrix = omatrix;
                    tempmatrix(endtick+1:end, YY ) = -60;
                    hsn = surf( XX, YY, tempmatrix( XX, YY)' );
                else
                    hsn = surf( XX(1:endtick), YY, omatrix( XX(1:endtick), YY)' );
                end
                axis([ 1 xmax+10 1 ymax+2 zmin zmax ]);
                set(gca,'Color', [0.8 0.8 0.8])
                set(gca,'tickdir', 'out')
                if k==1
                    colormap(cmap)
                    cm = colormap;
                    cm( end, : ) = [ 1 1 1 ];
                    colormap( cm );
                end
                view(viewazel)
                shading('interp')
                hsntext = sprintf('%s:  Vm Tracings, Cells %s', fname, cellrange );
                title( hsntext );
                xlabel('Time (ticks)')
                zlabel(['Hz, ' num2str(window/fsv) ' s window'] )
                M(k) = getframe(snfig);
            end
            %movie( M, 1, framerate );
            movie2avi( M, movfname, 'compression', 'indeo5', 'quality', 100, 'fps', framerate );
            fprintf(1,'\n--> Saved movie to file ''%s''.', movfname)
            %movieview(movfname)
            zoom on
        else % no movie
            hsn = surf( XX, YY, omatrix' );
            colormap(cmap)
            cm = colormap;
            cm( end, : ) = [ 1 1 1 ];
            colormap( cm );
            view(viewazel)
            shading('interp')
            hsntext = sprintf('%s:  Vm Tracings, Cells %s', fname, cellrange );
            title( hsntext );
            xlabel('Time (ticks)')
            zlabel(['Hz, ' num2str(window/fsv) ' s window'] )
            set(gca,'Color', [0.8 0.8 0.8])
            set(gca,'tickdir', 'out')
            zoom on
        end
    end
    
    
    if strcmp(plotmode,'raw2d')
        newpos = [ pos(1) pos(2)-pos(4)-75 pos(3) pos(4) ];
        cnfig = figure('units', 'normalized','position', figpos);
        if ~isempty(strmatch(play,{'movie', 'both'} ) )
            fprintf(1,'\n--> Making MOVIE of %d frames to play at %d /sec...', nframes, framerate)
            movfname = [ matname(1:end-4) '.avi' ];
            ticksperframe = ceil( size( omatrix, 1 )/nframes );
            for k = 1:nframes
                endtick = k*ticksperframe;
                if endtick > size( omatrix, 1 )
                    endtick = size( omatrix, 1 );
                end      
                if strcmp(basecolor,'on') % forces a -60 mV blue baseline rather than transparent
                    tempmatrix = omatrix;
                    tempmatrix(endtick+1:end, YY ) = -60;
                    hpn = pcolor( XX, YY, tempmatrix( XX, YY)'  );
                else
                    hpn = pcolor( XX(1:endtick), YY, omatrix( XX(1:endtick), YY)'  );
                end
                axis([ 1 xmax+10 1 ymax+2 ]);
                colormap( cmap )
                % adjust colormap to be shades of blue below threshold, white above
                rest = -65;
                minvolt = -80;
                maxvolt = spikepeak;
                cmaplow = round(minvolt):round(rest);
                cmapmid = round(rest)+1:round(thresh);
                cmaphi = round(thresh)+1:round(maxvolt);
                lenlow = length( cmaplow );
                lenmid = length( cmapmid );
                lenhi = length( cmaphi );
                lentotal = lenlow + lenmid + lenhi;
                bluevalues = linspace( .5, 1, length( cmaplow ) );
                greenvalues = linspace( 0, 1, length( cmapmid ) );
                whitevalues = zeros( 1, length(cmaphi) );
                zerocollow = zeros( lenlow, 1);
                zerocolmid = zeros( lenmid, 1);
                cnew = [ zerocollow, zerocollow, bluevalues' ];
                cnew = [ cnew; zerocolmid, greenvalues', ones( lenmid, 1) ];
                cnew = [ cnew; ones( lenhi, 3) ];
                colormap(cnew);    
                
                shading('flat')
                hpntext = sprintf('%s: Vm Tracings, Cells %s', fname, cellrange );
                title( hpntext )
                xlabel('Time (ticks)')
                
                % tick increment
                xt = get(gca,'xtick');
                currentincrem = xt(2) - xt(1);
                newincr = currentincrem/xtick2dfactor;
                xt = [xt(1):newincr:xt(end)];
                set(gca,'XTickMode','manual')
                set(gca, 'xtick', xt)
                
                set(gca,'Color', [0.8 0.8 0.8])
                set(gca,'tickdir', 'out')
                zoom on
                if k == nframes
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    figure(cnfig);
                    blim = [nanmin(omatrix(:)) nanmax(omatrix(:))];
                    
                    % colorbar
                    cmapLen = size(colormap,1);
                    hax_cbar = axes('pos',[.91 .275 .03 .525]);
                    himage_cbar = image([0 1],[0 1],(1:cmapLen)');
                    set(himage_cbar,'cdatamapping','scaled','erase','none');
                    set(hax_cbar, ...
                        'draw','fast', ...
                        'fontsize',8, ...
                        'box','on', ...
                        'xticklabel','', ...
                        'Ydir','normal', 'YAxisLocation','right', 'xtick',[]);
                    
                    % colorbar indicator
                    hax_cbar_ind = axes('pos',[.885+.01 .275 .015 .525]);
                    set(hax_cbar_ind,'vis','off','xlim',[0 1],'ylim',[0 1], ...
                        'draw','fast', ...
                        'fontsize',8, ...
                        'yaxisloc','right');
                    
                    % Update colorbar
                    set(himage_cbar, 'ydata',blim, 'cdata', (1:cmapLen)');
                    set(hax_cbar,'ylim',blim);
                    set(hax_cbar_ind, 'ylim',blim);
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                end
                drawnow
                M(k) = getframe(cnfig);
            end
            movie2avi( M, movfname, 'compression', 'indeo5', 'quality', 100, 'fps', framerate );
            fprintf(1,'\n--> Saved movie to file ''%s''.', movfname)
        else
            hpn = pcolor( XX, YY, omatrix' );
            colormap( cmap )
            % adjust colormap to be shades of blue below threshold, white above
            rest = -65;
            minvolt = nanmin(omatrix(:));
            maxvolt = spikepeak;
            cmaplow = round(minvolt):round(rest);
            cmapmid = round(rest)+1:round(thresh);
            cmaphi = round(thresh)+1:round(maxvolt);
            
            %debug
            lenlow = length( cmaplow );
            lenmid = length( cmapmid );
            lenhi = length( cmaphi );
            lentotal = lenlow + lenmid + lenhi;
            bluevalues = linspace( 0, 1, length( cmaplow ) );
            greenvalues = linspace( 0.25, .75, length( cmapmid ) );
            whitevalues = zeros( 1, length(cmaphi) );
            zerocollow = zeros( lenlow, 1);
            zerocolmid = zeros( lenmid, 1);
            cnew = [ zerocollow, zerocollow, bluevalues' ];
            cnew = [ cnew; zerocolmid, greenvalues', ones( lenmid, 1) ];
            cnew = [ cnew; ones( lenhi, 3) ];
            colormap(cnew);    
            
            shading('flat')
            hpntext = sprintf('%s: Vm Tracings, Cells %s', fname, cellrange );
            title( hpntext )
            xlabel('Time (ticks)')
            
            % tick increment
            xt = get(gca,'xtick');
            currentincrem = xt(2) - xt(1);
            newincr = currentincrem/xtick2dfactor;
            xt = [xt(1):newincr:xt(end)];
            set(gca,'XTickMode','manual')
            set(gca, 'xtick', xt)
            
            set(gca,'Color', [0.8 0.8 0.8])
            set(gca,'tickdir', 'out')
            zoom on
            
            if 1
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                figure(cnfig);
                blim = [nanmin(omatrix(:)) nanmax(omatrix(:))];
                % colorbar
                cmapLen = size(colormap,1);
                hax_cbar = axes('pos',[.91 .275 .03 .525]);
                himage_cbar = image([0 1],[0 1],(1:cmapLen)');
                set(himage_cbar,'cdatamapping','scaled','erase','none');
                set(hax_cbar, ...
                    'draw','fast', ...
                    'fontsize',8, ...
                    'box','on', ...
                    'xticklabel','', ...
                    'Ydir','normal', 'YAxisLocation','right', 'xtick',[]);
                
                % colorbar indicator
                hax_cbar_ind = axes('pos',[.885+.01 .275 .015 .525]);
                set(hax_cbar_ind,'vis','off','xlim',[0 1],'ylim',[0 1], ...
                    'draw','fast', ...
                    'fontsize',8, ...
                    'yaxisloc','right');
                
                % Update colorbar
                
                set(himage_cbar, 'ydata', blim, 'cdata', (1:cmapLen)');
                set(hax_cbar,'ylim', blim);
                set(hax_cbar_ind, 'ylim', blim);
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            end
            drawnow
        end
    end
    if ~isempty(strmatch(play,{'sound', 'both'} ) )
        soundon( matname, omatrix, fsv, soundfreq, soundsubsample ); 
    end
end

% Vm TRACING SUBPLOTS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(plotmode,'rawsubplots')
    spfig = figure;
    nplots = size( omatrix, 2 );
    sptext = sprintf('%s: Vm Subplots, Cells %s', fname, cellrange );
    for p = 1:nplots
        subplot( nplots, 1, p )
        if p == 1, title(sptext), end; hold on
        plot( omatrix( :, nplots-p+1), '-', 'color', rand(1,3) )
        set(gca,'tickdir', 'out')
        hold on
        zoom on
        % tick increment
        xt = get(gca,'xtick');
        currentincrem = xt(2) - xt(1);
        newincr = currentincrem/xtick2dfactor;
        xt = [xt(1):newincr:xt(end)];
        set(gca,'XTickMode','manual')
        set(gca, 'xtick', xt)
        
        %grid on
        %ylabel( collabels(nplots-p+1) )
    end
    if ~isempty(strmatch(play,{'sound', 'both'} ) )
        soundon( matname, omatrix, fsv, soundfreq, soundsubsample ); 
    end
end

% DOT RASTER PLOT  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(plotmode,'rawdot')
    spikeindices = find( omatrix == spikepeak);
    if ~isempty( spikeindices )
        dfig = figure;
        I = []; J = [];
        fprintf(1,'\n--> EXTRACTING points for dot raster...')
        for dot = 1:length(spikeindices)
            [ i, j ] = ind2sub( size(omatrix), spikeindices( dot ));
            I = [ I, i ]; J = [ J, j ];
        end
        xmax = size( omatrix, 1 ); ymax = size( omatrix, 2 );
        
        if ~isempty(strmatch(play,{'movie', 'both'} ) )
            fprintf(1,'\n--> Making MOVIE of %d frames to play at %d /sec...', nframes, framerate)
            movfname = [ matname(1:end-4) '.avi' ];
            ticksperframe = ceil( size( omatrix, 1 )/nframes );
            for k = 1:nframes
                endtick = k*ticksperframe;
                if endtick > size( omatrix, 1 )
                    endtick = size( omatrix, 1 );
                end                    
                xindices = find( I <= endtick );
                hp = plot( I(xindices), J(xindices), 'k.' );
                set(hp, 'MarkerSize', dotsize )
                axis([ 1 xmax+10 1 ymax+2 ]);
                % tick increment
                xt = get(gca,'xtick');
                currentincrem = xt(2) - xt(1);
                newincr = currentincrem/xtick2dfactor;
                xt = [xt(1):newincr:xt(end)];
                set(gca,'XTickMode','manual')
                set(gca, 'xtick', xt)
                
                M(k) = getframe;
            end
            movie2avi( M, movfname, 'compression', 'indeo5', 'quality', 100, 'fps', framerate );
            fprintf(1,'\n--> Saved movie to file ''%s''.', movfname)
            %movie( M, 1, framerate );
            %movieview(movfname)
            zoom on
        else % no movie
            hp = plot( I, J, 'k.' );
            set(hp, 'MarkerSize', dotsize )
            axis([ 1 xmax+10 1 ymax+2 ]);
            % tick increment
            xt = get(gca,'xtick');
            currentincrem = xt(2) - xt(1);
            newincr = currentincrem/xtick2dfactor;
            xt = [xt(1):newincr:xt(end)];
            set(gca,'XTickMode','manual')
            set(gca, 'xtick', xt)
            
            zoom on
        end
        set(gca,'tickdir', 'out')
        if ~isempty(strmatch(play,{'sound', 'both'} ) )
            smatrix = sparse( zeros( size( omatrix ) ) );
            smatrix( spikeindices ) = 1;
            soundon( matname, smatrix, fsv, soundfreq, soundsubsample ); 
        end           
    else
        fprintf(1,'\n--> NO SPIKES FOUND.')
    end
    dtext = sprintf('%s: Spike Dot Raster, Cells %s', fname, cellrange );
    title( dtext )
    xlabel('Time (ticks)')
    ylabel('Cell ID')
    zoom on
end

% DENSITY PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(plotmode,'density')
    density( fname, omatrix, tickcompress, cellcompress, spikepeak, fsv, xtick2dfactor,figpos );
end
% AVG PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isempty(strmatch(plotmode,{'avg2d', 'avg3d'} ) )
    matrix = window_it( omatrix, window,...
        refrange, comprange, spikepeak, fsv );
    
    %debug
    size(matrix)
    mean(matrix(:))

    % sort cells by change in rate from ref to comparison ranges
    if ~strcmp( comprange, 'nosort')
        [ matrix, I ] = sortcomp( matrix, omatrix, refrange, comprange );
        celllist = [1:size(omatrix, 2)];
        %order = eval( [ 'celllist(' cellrange ')' ] ); % numeric
        order = celllist( I ); 
    end
    matrix = [zeros( 1, size( matrix, 2)); matrix ]; % add a row
    matrix = [matrix  zeros( size( matrix, 1), 1)]; % add a col
    pos = [155         690        1395         420];
    XX = [1:size( matrix, 1)];
    YY = [1:size( matrix, 2)];
    xmax = size( matrix, 1) ; ymax = size( matrix, 2 ); 
    zmin = nanmin( matrix(:)); zmax = nanmax( matrix(:));
    if strcmp(plotmode,'avg3d')
        sfig = figure('units', 'normalized','position', figpos);        
        if ~isempty(strmatch(play,{'movie', 'both'} ) )
            fprintf(1,'\n--> Making MOVIE of %d frames to play at %d /sec...', nframes, framerate)
            movfname = [ matname(1:end-4) '.avi' ];
            ticksperframe = ceil( size( matrix, 1 )/nframes );
            for k = 1:nframes
                endtick = k*ticksperframe;
                if endtick > size( matrix, 1 )
                    endtick = size( matrix, 1 );
                end       
                
                if strcmp(basecolor,'on') % forces a -60 mV blue baseline rather than transparent
                    tempmatrix = matrix;
                    tempmatrix(endtick+1:end, YY ) = 0;
                    hs = surf( XX, YY, tempmatrix( XX, YY)' );
                else
                    hs = surf( XX(1:endtick), YY, matrix( XX(1:endtick), YY)' );
                end
                hold on
                %figure(sfig)
                % plot( [ 11000 11000], [ 0 80], 'c-')
                axis([ 1 xmax+10 1 ymax+2 zmin  zmax ]);
                if k==1
                    colormap( cmap )
                    cm = colormap;
                    cm( end, : ) = [ 1 1 1 ];
                    colormap( cm );
                end
                view(viewazel)
                % white for nanmax of colormap
                shading('interp')
                stext = sprintf('%s: %d-Tick Windowed Rates, Cells %s', fname, window, cellrange );
                title( stext )
                xlabel('Time (ticks)')
                zlabel(['Hz, ' num2str(window/fsv) ' s window'] )
                % resort cell id labels also
                if ~strcmp( refrange, 'nosort')
                    set(gca,'YTickLabel', order ) 
                end
                set(gca,'Color', [0.8 0.8 0.8])
                M(k) = getframe(sfig);
            end
            %movie( M, 1, framerate );
            movie2avi( M, movfname, 'compression', 'indeo5', 'quality', 100, 'fps', framerate );
            fprintf(1,'\n--> Saved movie to file ''%s''.', movfname)
            %movieview(movfname)
            zoom on
        else % no movie
            hs = surf( XX, YY, matrix' );
            colormap( cmap )
            cm = colormap;
            cm( end, : ) = [ 1 1 1 ];
            colormap( cm );
            view(viewazel)
            % white for nanmax of colormap
            shading('interp')
            stext = sprintf('%s: %d-Tick Windowed Rates, Cells %s', fname, window, cellrange );
            title( stext )
            xlabel('Time (ticks)')
            zlabel(['Hz, ' num2str(window/fsv) ' s window'] )
            % resort cell id labels also
            if ~strcmp( refrange, 'nosort')
                set(gca,'YTickLabel', order ) 
            end
            set(gca,'Color', [0.8 0.8 0.8])
            set(gca,'tickdir', 'out')
            zoom on
        end
    end
    
    if strcmp(plotmode,'avg2d')
        newpos = [ pos(1) pos(2)-pos(4)-75 pos(3) pos(4) ];
        cfig = figure('units', 'normalized', 'position', figpos);
        if ~isempty(strmatch(play,{'movie', 'both'} ) )
            fprintf(1,'\n--> Making MOVIE of %d frames to play at %d /sec...', nframes, framerate)
            movfname = [ matname(1:end-4) '.avi' ];
            ticksperframe = ceil( size( matrix, 1 )/nframes );
            for k = 1:nframes
                endtick = k*ticksperframe;
                if endtick > size( matrix, 1 )
                    endtick = size( matrix, 1 );
                end                    
                hp = pcolor( XX(1:endtick), YY, matrix( XX(1:endtick), YY)' );
                axis([ 1 xmax+10 1 ymax+2 ]);
                colormap( cmap )
                cm = colormap;
                cm( end, : ) = [ 1 1 1 ];
                colormap( cm );
                shading('flat')
                ctext = sprintf('%s: %d-Tick Windowed Rates, Cells %s', fname, window, cellrange );
                title( ctext )
                xlabel('Time (ticks)')
                % tick increment
                xt = get(gca,'xtick');
                currentincrem = xt(2) - xt(1);
                newincr = currentincrem/xtick2dfactor;
                xt = [xt(1):newincr:xt(end)];
                set(gca,'XTickMode','manual')
                set(gca, 'xtick', xt)
                
                if ~strcmp( refrange, 'nosort')
                    set(gca,'YTickLabel', order ) 
                end
                set(gca,'Color', [0.8 0.8 0.8])
                set(gca,'tickdir', 'out')
                zoom on
                if k == nframes  % add colorbar legend
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    figure(cfig);
                    blim = [nanmin(matrix(:)) nanmax(matrix(:))];
                    
                    % colorbar
                    cmapLen = size(colormap,1);
                    hax_cbar = axes('pos',[.91 .275 .03 .525]);
                    himage_cbar = image([0 1],[0 1],(1:cmapLen)');
                    set(himage_cbar,'cdatamapping','scaled','erase','none');
                    set(hax_cbar, ...
                        'draw','fast', ...
                        'fontsize',8, ...
                        'box','on', ...
                        'xticklabel','', ...
                        'Ydir','normal', 'YAxisLocation','right', 'xtick',[]);
                    
                    % colorbar indicator
                    hax_cbar_ind = axes('pos',[.885+.01 .275 .015 .525]);
                    set(hax_cbar_ind,'vis','off','xlim',[0 1],'ylim',[0 1], ...
                        'draw','fast', ...
                        'fontsize',8, ...
                        'yaxisloc','right');
                    
                    % Update colorbar
                    set(himage_cbar, 'ydata',blim, 'cdata', (1:cmapLen)');
                    set(hax_cbar,'ylim',blim);
                    set(hax_cbar_ind, 'ylim',blim);
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                end
                M(k) = getframe(cfig);
            end
            %movie( M, 1, framerate );
            movie2avi( M, movfname, 'compression', 'indeo5', 'quality', 100, 'fps', framerate );
            fprintf(1,'\n--> Saved movie to file ''%s''.', movfname)
            %movieview(movfname)
            zoom on
            
        else % no movie
            
            hp = pcolor( XX, YY, matrix' );
            colormap( cmap )
            
            cm = colormap;
            cm( end, : ) = [ 1 1 1 ];
            colormap( cm );
            shading('flat')
            ctext = sprintf('%s: %d-Tick Windowed Rates, Cells %s', fname, window, cellrange );
            title( ctext )
            xlabel('Time (ticks)')
            % tick increment
            xt = get(gca,'xtick');
            currentincrem = xt(2) - xt(1);
            newincr = currentincrem/xtick2dfactor;
            xt = [xt(1):newincr:xt(end)];
            set(gca,'XTickMode','manual')
            set(gca, 'xtick', xt)
            
            if ~strcmp( refrange, 'nosort')
                set(gca,'YTickLabel', order ) 
            end
            set(gca,'Color', [0.8 0.8 0.8])
            set(gca,'tickdir', 'out')
            zoom on
            if 1  % add colorbar legend
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                figure(cfig);
                blim = [nanmin(matrix(:)) nanmax(matrix(:))];
                
                % colorbar
                cmapLen = size(colormap,1);
                hax_cbar = axes('pos',[.91 .275 .03 .525]);
                himage_cbar = image([0 1],[0 1],(1:cmapLen)');
                set(himage_cbar,'cdatamapping','scaled','erase','none');
                set(hax_cbar, ...
                    'draw','fast', ...
                    'fontsize',8, ...
                    'box','on', ...
                    'xticklabel','', ...
                    'Ydir','normal', 'YAxisLocation','right', 'xtick',[]);
                
                % colorbar indicator
                hax_cbar_ind = axes('pos',[.885+.01 .275 .015 .525]);
                set(hax_cbar_ind,'vis','off','xlim',[0 1],'ylim',[0 1], ...
                    'draw','fast', ...
                    'fontsize',8, ...
                    'yaxisloc','right');
                
                % Update colorbar
                set(himage_cbar, 'ydata',blim, 'cdata', (1:cmapLen)');           
                set(hax_cbar,'ylim',blim);
                set(hax_cbar_ind, 'ylim',blim);
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            end
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get_list_of_fnames()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function flist = get_list_of_fnames(txt_fname, starindex)

d = dir('.');
d(1:2)=[];
ld = length(d);
flist = {};
lt = length(txt_fname);
firstpart = txt_fname(1:(starindex-1));
lastpart = txt_fname((starindex+1):end);
for f = 1:ld
    %d(f).name
    if findstr( d(f).name, firstpart ) & findstr( d(f).name, lastpart )
        flist{length(flist)+1} = d(f).name;
    end
end
flist = sort( flist );
if isempty(flist)
    emsg = sprintf('!!! No files matching "%s*" were found', firstpart );
    error( emsg )
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% test_tick_col()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function matrix = test_tick_col( matrix, celldata)

% if cell data in rows, first rotate the matrix
if celldata(1) == 'r'
    matrix = matrix';
end
% detect and skip first column (row, if celldata 'row')  if it is an integer label for the tick number
if matrix( 2, 1 ) > matrix( 1, 1) & matrix( 3, 1 ) > matrix( 2, 1)
    if round( matrix( 1, 1 ) ) == matrix( 1, 1 ) & round( matrix( 2, 1 ) ) == matrix( 2, 1 ) & round( matrix( 3, 1 ) ) == matrix( 3, 1 )
        matrix = matrix(:,2:end);
    end
end

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% window_it()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function winavg = window_it( matrix, tickwindow,...
    refrange, comprange, spikepeak, fsv );

ncells =  size(matrix, 2);
maxtick = size(matrix, 1);
nwindows = maxtick-tickwindow;
zmatrix = repmat( 0, size(matrix) );
zmatrix( find( matrix == spikepeak ) ) = 1;
clear matrix
fprintf(1,'\n--> %% done:')
total_cellticks = ncells*nwindows;
fivepctdone = round(0.05*total_cellticks);
winavg = repmat(nan, nwindows, ncells); % initialize for efficiency
if 0
    for cell = 1:ncells
        for tick = 1:nwindows
            thiswindow = zmatrix( tick:tick+tickwindow, cell );
            winavg( tick, cell ) = sum( thiswindow );
            
            if ~rem( (cell-1)*nwindows + tick, fivepctdone )
                fprintf(1,' %0.0f ', round(100*((cell-1)*nwindows + tick)/total_cellticks))
            end
        end
    end
else
    for cell = 1:ncells
        winavg( 1, cell ) = sum( zmatrix( 1:tickwindow, cell ) );
        
        for tick = 2:nwindows
            prevendvalue = zmatrix( (tick-1), cell );
            currentendvalue = zmatrix( tick+tickwindow, cell );
            
            if prevendvalue | currentendvalue
                if prevendvalue & currentendvalue
                    winavg( tick, cell ) = winavg( tick-1, cell );
                elseif prevendvalue
                    addto = -1;
                    winavg( tick, cell ) = winavg( tick-1, cell )-1;
                else % must be only new spike entering, none leaving window
                    winavg( tick, cell ) = winavg( tick-1, cell )+1;
                end
            else 
                winavg( tick, cell ) = winavg( tick-1, cell );
            end
            lastwindow = zmatrix( tick:tick+tickwindow, cell );
            if ~rem( (cell-1)*nwindows + tick, fivepctdone )
                fprintf(1,' %0.0f ', round(100*((cell-1)*nwindows + tick)/total_cellticks))
            end
        end
    end
    
end
winavg = winavg * fsv/tickwindow; % convert to Hz


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sortcomp()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ matrix, I ]  = sortcomp( matrix, omatrix, refrange, comprange )

ncells =  size(matrix, 2);
for cell = 1:ncells
    % range
    comprate = eval( [ 'mean( omatrix(' comprange ', cell) );' ]);
    if strcmp( refrange, 'nosort') % absolute values
        ratechange(cell) = comprate;
    else  % compare comprange to a refrange
        refrate = eval( [ 'mean( omatrix(' refrange ', cell) );' ]);
        ratechange(cell) = comprate - refrate;
    end
end
[ Y, I ] = sort(ratechange);

matrix = matrix( :, I );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parsearg()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ tokenlist, valuelist ] = parsearg( args, legal_tokens )
% examine each arg, test if a token; if so, return value; else return as an arg

nargs = length( args );
if rem(nargs,2) ~= 0 % not even number
    error('Odd number of arguments: must be in pairs as (...,''parameter'', ''value'',...)')
    return
end

vindex = 0;
for a = 1:2:nargs
    vindex = vindex + 1;
    if ~any( strcmp( args{a}, legal_tokens ) )
        error([ args{a} ' is not a legal parameter.' ])
    else
        tokenlist{vindex} = args{a};
        valuelist{vindex} = args{a+1};
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% density()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Density = density( fname, submatrix, colcompression, rowcompression, spike_peak_mV, fsv,xtick2dfactor,figpos )
% creates density plot for spikes   rev 15 Jly 01  PHG
%
% EX: density( 'oneyes.mat', 90, 90, 40 );
%
% matfile:      .mat filename that stores the object called submatrix (this is saved by plotmany.m)
% colcompression:  compression factor for columns (ticks); recommend 50-100
% rowcompression:  compression factor for columns (ticks); recommend 50-100
% spike_peak_mV:   used to count number of spikes in a region 

submatrix = submatrix';
if nanmax( submatrix(:) ) ~= 1  % not binarized yet
    ind = find( submatrix == spike_peak_mV );
    submatrix = zeros( size( submatrix) );
    submatrix(ind ) = 1;
end

n_rowbins = floor( size(submatrix,1)/rowcompression);
n_colbins = floor( size(submatrix,2)/colcompression);

if n_rowbins < 2
    error('Less than 2 cell bins: make ''cellcompress'' smaller [DEFAULT is 50]')
end
if n_colbins < 2
    error('Less than 2 tick bins: make ''tickcompress'' smaller [DEFAULT is 50]')
end
Density = [];
for rbin = 1:n_rowbins
    for cbin= 1:n_colbins   
        %[rbin,cbin]
        subregion = submatrix( ((rbin-1)*rowcompression+1):rbin*rowcompression, ((cbin-1)*colcompression+1):cbin*colcompression );
        
        Density( rbin, cbin) = sum( subregion(:) );
    end
end

% rescale Density to reflect spikes/cell/sec
Density = Density / (colcompression/fsv) / rowcompression;
densfig =figure;
[ matrix, hc ]  = contourf(Density,30);
set(gcf,'units', 'normalized','position', figpos );

% tick increment
xt = get(gca,'xtick');
currentincrem = xt(2) - xt(1);
newincr = currentincrem/xtick2dfactor;
xt = [xt(1):newincr:xt(end)];
set(gca,'XTickMode','manual')
set(gca, 'xtick', xt)

xtext = [ 'Time Bin (each ' num2str(colcompression) ' ticks)'];
ytext = [ 'Cell Bin (each ' num2str(rowcompression) ' cells)'];
shading('flat')

title([fname ': DENSITY OF FIRING RATES (mean per cell per second)'])
xlabel(xtext)
ylabel(ytext)
blim = [nanmin(Density(:)) nanmax(Density(:))];
set(gca,'tickdir', 'out')

% colorbar
cmapLen = size(colormap,1);
hax_cbar = axes('pos',[.91 .275 .03 .525]);
himage_cbar = image([0 1],[0 1],(1:cmapLen)');
set(himage_cbar,'cdatamapping','scaled','erase','none');
set(hax_cbar, ...
    'draw','fast', ...
    'fontsize',8, ...
    'box','on', ...
    'xticklabel','', ...
    'Ydir','normal', 'YAxisLocation','right', 'xtick',[]);

% colorbar indicator
hax_cbar_ind = axes('pos',[.885+.01 .275 .015 .525]);
set(hax_cbar_ind,'vis','off','xlim',[0 1],'ylim',[0 1], ...
    'draw','fast', ...
    'fontsize',8, ...
    'yaxisloc','right');

% Update colorbar
%XXX UINT8 change:
set(himage_cbar, 'ydata',blim, 'cdata', (1:cmapLen)');
set(hax_cbar,'ylim',blim);
set(hax_cbar_ind, 'ylim',blim);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% soundmovie()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function soundon( matname, matrix, fsv, soundfreq, soundsubsample )

if issparse(matrix), matrix = full( matrix ); end
wavfname = [ matname(1:end-4) '.wav' ];

if strcmp( soundsubsample, 'random' )
    % zero out all but one cell each tick
    ncells = size( matrix, 2 );
    Cells = 1:ncells;
    for tick = 1:size( matrix, 1 )
        randcell = ceil( ncells*rand );
        matrix( tick, setdiff( Cells, randcell ) ) = 0;
    end
else
    matrix = eval( [ 'matrix(:,' soundsubsample ');'] );% subsample of cells
end

if size(matrix, 2) > 1
    smatrix= nanmax(matrix'); % compress all spikes into one line
    wavplay(smatrix, soundfreq,'sync')
    wavwrite( smatrix, soundfreq, wavfname );
    clear smatrix
else
    wavplay(matrix,soundfreq,'sync')
    wavwrite( matrix, soundfreq, wavfname );
end
fprintf(1,'\n->Saved SOUND file ''%s''', wavfname);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% checkfornans()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fdata = checkfornans( fname, fdata, deletenan )

[I,J] = ind2sub( size(fdata), find( isnan(fdata) ) );
rows_to_delete = [];
if ~isempty( I )  % at least one nan found
    fprintf(1,'\n!!! FYI, %d NANs found in ''%s'' at:', length(I), fname);
    cellold = 0;
    nticksinarow = 0;
    for nansub = 1:length(I)
        if J(nansub) ~= cellold 
            if deletenan(1) == 'y'
                rows_to_delete = [ rows_to_delete J(nansub) ];
            end                    
            if nticksinarow > 1
                fprintf(1,'\n\t\t[%d ticks in row for cell above]', nticksinarow );
            end
            fprintf(1,'\n\t(  cell %d , tick %d)',  J(nansub), I( nansub) );
            cellold = J(nansub);
            nticksinarow = 1;
        else
            nticksinarow = nticksinarow + 1;
        end
    end
    uniquecellids =unique(J);
    if J(nansub) == uniquecellids(end);
        fprintf(1,'\n\t\t[%d ticks in row for cell above]', nticksinarow );
    end
    if deletenan(1) == 'y'
        fdata( :, rows_to_delete ) = [];
        fprintf(1,'\n--> Deleted %d rows containing NaNs.', length( rows_to_delete ) );
    end
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% slashunderscore()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fname= slashunderscore( fname )
% for titles, escape '_' titles to prevent underscore

slashindices = find( fname == '_' );

if ~isempty( slashindices )
    for sindx = 1:length(slashindices)
        nextindex = find(fname == '_' );
        firstpart = fname(1:nextindex(sindx)-1);
        endpart = fname(nextindex(sindx):end);
        fname = [ firstpart '\' endpart ];
    end
else
    fname=fname;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cf()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cf( )
% close all figures

figure(1); %select figure 1
while (gpf ~=1) %close all
    close(gpf);
end
close(1)
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  END OF LEARNPLOT.M  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%