% mytest.m
% Calls MultiClassifyByDiffusion (engine), Diffuser_To_Neighbors (diffusion), Classifier_MultiClass_PL (classifier)
% Ann Lee, March 2004

% DATA SET
set=3; %3
switch set
    case 1 % annulus + center
        load('workspace.mat');
        cPoints=cData; clear v*
        N=size(cPoints,1); figure(10), 
        P=spalloc(2,N,2*N); P(1,1)=1; P(2,2000)=1; NumClasses=size(P,1),     % P(NumClasses,N)
        Nc = full(sum(P>0,2))';  % number of labeled points in each class
        cDiffuserExtraParam.Heat=P;
        cDiffuserExtraParam.Delta=0.4; 
        cDiffuserExtraParam.Precision=0.1;
        figure(1), clf, scatter(cPoints(:,1),cPoints(:,2),15,'y.'); drawnow,hold  on; 
        tmpidx=find(sum(P,1)>0); [dummy,tmplabels]=max(P,[],1); scatter(cPoints(tmpidx,1),cPoints(tmpidx,2),15,tmplabels(tmpidx)); hold on, pause,
    case 2 % two crossing curves
        cPoints = generate2d('cubic+line',0,1000); N=size(cPoints,1); tmpidx=[1 1000 1001 2000]; 
        P=spalloc(2,N,2*N); P(1,1)=1/2; P(1,1000)=1/2;
        P(2,1001)=1/2; P(2,2000)=1/2; 
        NumClasses=size(P,1);     % P(NumClasses,N)
        Nc = full(sum(P>0,2))';  % number of labeled points in each class
        %clf, figure(10), plot(cPoints(:,1), cPoints(:,2), 'y.', cPoints(tmpidx,1), cPoints(tmpidx,2),'rx'); 
        %axis([-1 1 -1 1]); drawnow; %hold on;
        cDiffuserExtraParam.Heat=P;
        cDiffuserExtraParam.Delta=0.04; %1; 
        cDiffuserExtraParam.Precision=0.1;
        figure(10), clf, scatter(cPoints(:,1),cPoints(:,2),5,'y'); drawnow, hold  on;
        tmpidx=find(sum(P,1)>0); [dummy,tmplabels]=max(P,[],1); scatter(cPoints(tmpidx,1),cPoints(tmpidx,2),15,tmplabels(tmpidx)); hold on, pause,
    case 3 % mixed set with 5 clusters
        load mix_data.mat;   % 5 clusters
        cPoints=x; N=size(cPoints,1); 
        axis([-2 7 -2 7]); drawnow; %%hold  on;
        NumClasses=length(idx); tmpind=sub2ind([NumClasses,N],1:NumClasses,idx);
        P=spalloc(NumClasses,N,NumClasses*N); P(tmpind)=1;    % P(NumClasses,N)
        Nc = full(sum(P>0,2))';  % number of labeled points in each class
        cDiffuserExtraParam.Heat=P;
        cDiffuserExtraParam.Delta=0.4; 
        cDiffuserExtraParam.Precision=0.1;
        figure(1), clf, scatter(cPoints(:,1),cPoints(:,2),5,'y'); hold  on;
        tmpidx=find(sum(P,1)>0); [dummy,tmplabels]=max(P,[],1); scatter(cPoints(tmpidx,1),cPoints(tmpidx,2),40,tmplabels(tmpidx)); pause,
    case 4 % corona data
        %pNumberOfPoints= 500;
        % Generate the points
        %cPoints=GenerateSampleCoronaData(pNumberOfPoints,1,4,0.4);
        %cDiffuserExtraParam.Heat=spalloc(2,pNumberOfPoints,2*pNumberOfPoints);
        %cDiffuserExtraParam.Heat(1,1)=1; cDiffuserExtraParam.Heat(2,500)=1;
        %cDiffuserExtraParam.Delta=0.5;
        %cDiffuserExtraParam.Precision=0.1;
        N=2000; NumClasses=2;
        cPoints=GenerateSampleCoronaData(N,1,5,0.6); %0.8);
        %cPoints=GenerateSampleCoronaData(N,1,5,0.2); %0.8);
        cDiffuserExtraParam.Heat=spalloc(NumClasses,N,NumClasses*N);  % NEW: NumClasses X N
        cDiffuserExtraParam.Heat(1,1)=1; cDiffuserExtraParam.Heat(2,N)=1; % Note: one labeled point from each class
        cDiffuserExtraParam.Delta=0.2; %0.6
        cDiffuserExtraParam.Precision=0.1;
        figure(1), scatter(cPoints(:,1),cPoints(:,2),5,'y.'); hold  on; P=cDiffuserExtraParam.Heat;
        tmpidx=find(sum(P,1)>0); [dummy,tmplabels]=max(P,[],1); scatter(cPoints(tmpidx,1),cPoints(tmpidx,2),15,tmplabels(tmpidx)); hold on, pause,
    case 5
        cPoints=(1:10)';
        cDiffuserExtraParam.Heat=spalloc(1,10,10); cDiffuserExtraParam.Heat(1)=1;
        cDiffuserExtraParam.Delta=1; cDiffuserExtraParam.Precision=0.01;
end

%------------------------------------------------------------------
% NN NEIGHBOR CLASSIFICATION (FOR COMPARISON)
%------------------------------------------------------------------
if 0
    disp('Nearest Neighbor Classification');
    tmpidx=find(sum(P,1)>0); [dummy,tmplabels]=max(P,[],1); 
    %scatter(cPoints(tmpidx,1),cPoints(tmpidx,2),40,tmplabels(tmpidx)); break,
    [dummy,tmpidx2]=sort(tmplabels(tmpidx));
    dtmp = dist2(cPoints(tmpidx(tmpidx2),:),cPoints);
    [tmp,tmpL] = min(dtmp,[],1);
    figure, scatter(cPoints(:,1),cPoints(:,2),15,tmpL);
    axis image, colorbar, title('1-NN Classification');
end


% Classifier parameters (Input to Classifier_BdryGrowth)
cClassifierExtraParam.Type='standard';
cClassifierExtraParam,

global gClassifierObj 
global gDiffuseObj

% New parameters to "MultiClassifyByDiffusion", AL 3/16/05
cDiffuserExtraParam.GetWeightsFcn = @GetWeights_IsotropicGaussian;  
cDiffuserExtraParam.maxtime=201; GFlag=0; %0; 
cDiffuserExtraParam.Metric=cell(N,1); % new
cDiffuserExtraParam.MetricFlags = uint8(zeros(N,1));
cDiffuserExtraParam.Epsilon=2*cDiffuserExtraParam.Delta; %new
cDiffuserExtraParam.NormalizationType='Beltrami'; % NEW
cDiffuserExtraParam.NNsymm='ave'; % NEW
if GFlag==1  % add extra parameters for global diffusion
    cDiffuserExtraParam.EigFlag=1;  % cExtraParm.EigFlag=1;
    cDiffuserExtraParam.kNN=50;  % 
    cDiffuserExtraParam.MaxEigenVals=50; %100 %50;
    cDiffuserExtraParam.EigEnergyThres=0.99;
    if cDiffuserExtraParam.EigFlag==1 % can skip time steps
        disp('Global diffusion (with eigenvectors)')
        cDiffuserExtraParam.tstep=10; %10;
    else % step-wise global diffusion
        cDiffuserExtraParam.tstep=1;
        disp('Global diffusion (step-wise)')
    end
    % Call diffusion routine!
    MultiClassifyByDiffusion(cPoints,[],@Diffuser_GlobalEig,cDiffuserExtraParam,@Classifier_MultiClass_PL,cClassifierExtraParam);
else  % local diffusion (step-wise)
    disp('Local diffusion');
    cDiffuserExtraParam.tstep=1;
    % Call diffusion routine!
    MultiClassifyByDiffusion(cPoints,[],@Diffuser_ToNeighbors,cDiffuserExtraParam,@Classifier_MultiClass_PL,cClassifierExtraParam);
end
cDiffuserExtraParam,

% Do the diffusion
%ClassifyByDiffusion(cPoints,[],@Diffuser_ToNeighbors,cDiffuserExtraParam,@Classifier_BdryGrowth,cClassifierExtraParam ) 
%MultiClassifyByDiffusion(cPoints,[],@Diffuser_ToNeighbors,cDiffuserExtraParam,@Classifier_MultiClass_PL,cClassifierExtraParam);
%MultiClassifyByDiffusion(cPoints,[],@Diffuser_GlobalEig,cDiffuserExtraParam,@Classifier_MultiClass_PL,cClassifierExtraParam);
