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 MaxSpecTable MaxSpeedTable
close all
%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 = 23;

ngrp = 4;
grp{1} = 1:5  ; 
grp{2} = 6:13 ;
grp{3} = 14:19 ;
grp{4} = 20:23 ;

doHull=1; %if 1 does convex hull plot
doANOVA=1;

%A = squeeze(dVictor(5,:,:));
A = (1 - MaxSpeedTable).^1 ;
anovadim = 1;
plotdim1 	= 1;		% which dimension on x-axis
plotdim2 	= 3;		% 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;
A
% ---------------------------------------------------------------------------------
%    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 );
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 ));
   
   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

if ~doANOVA
    return
end

%anova stuff
%build a 9 x 4 matrix of Config(i,1)
ainput = nan * ones(9,4);

for i=1:ngrp
        pts = grp{i};
        for j =1:length(pts)
            ainput(j,i) = Config(pts(j),anovadim);
        end
    end 
    
    [p,table,stats] = anova1(ainput);
    [c,m] = multcompare(stats,.05,'off','bonferroni')


%===STOP HERE===============
return
%===========================

%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')