domdsprofile.m

MULTIDIMENSIONAL SCALING in matlab by Mark Steyvers 1999 ％needs optimization toolbox ％Modified b

function domds 
% ----------------------------------------------------------------------------------
%    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
--
% ----------------------------------------------------------------------------------
global dVictor ngrp grp dRossum tSpike
%no 			= 10;						% size of input matrix
%A  			= rand( no , no );	% dissimilarities

% no = 3;
% A = [ 0 3 4 ; ...
%       3 0 5; ...
%       4 5 0 ];

no = 20;
%A = squeeze(dVictor(5,:,:));
A = squeeze(dRossum(4,:,:));

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 );
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 );
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*' );
   grid on;
   axis equal;
   xlabel( sprintf( 'dimension %d' , plotdim1 ));
   ylabel( sprintf( 'dimension %d' , plotdim2 ));
   
   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 = vga; %jet(ngrp);
    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

%get the temporal profiles

%bin each spike train to form R, then append a col of ones 
%(for translation)
%time 50 to 100 mSec: centers 55 65 75 85 95
nbins = 5;
edgebins = [.25 .32 .4 .6 .8 ];
centerbins = [55 65 75 85 95]*.001;
for i=1:no
    R(i,1:nbins) = hist(tSpike{i},nbins);
   %R(i,1:nbins) = histc(tSpike{i},edgebins);
end

R(:,nbins+1) = 1;

% Config array from above is C
%compute P  such that C=RP
P = pinv(R)*Config ;
%delete the last row (cood of null train)
P = P(1:end-1,:);

%Plot profiles
figure(6)
clf
plot(P)
legend('d1','d2','d3','d4','d5','d6')
xlabel('Bin number')
ylabel('Weight in dimension d')