% Test klinreg on large data sets
%
% You must have klinreg_matlab installed and in the PATH
% to run this test.
%
clear all;
close all

DO_OTHERS = 0; 	% if 0, use only klinreg, otherwise compare with other methods
		% this requires the code for the other methods (not included in this package)


N = 10000;	% number of data

Npar=3;		% dimension p
Nmodes=2;	% number of modes 
Nexpe=100;	% number of repeated experiments

sigmaNoise = 0.2;	% noise standard deviation

Wexpe=[];
Wtmp = zeros(Npar,Nmodes);
Wtmp_ordered = zeros(Npar,Nmodes);
Y = zeros(N,1);
NMSEparam = zeros(Nexpe,1);
ClassifErrRate = zeros(Nexpe,1);
comp_time = zeros(Nexpe,1);

if DO_OTHERS
	NMSEparamME = zeros(Nexpe,1);
	ClassifErrRateME = zeros(Nexpe,1);
	comp_timeME = zeros(Nexpe,1);

	NMSEparamPE = zeros(Nexpe,1);
	ClassifErrRatePE = zeros(Nexpe,1);
	comp_timePE = zeros(Nexpe,1);
end

for expe=1:Nexpe
    expe
    
    % data
    X = rand(N,Npar)*10-5;
    Wt=zeros(Npar,Nmodes);
    while(sum((Wt(:,1)-Wt(:,2)).^2)<1)
        if (sum((Wt(:,1)-Wt(:,2)).^2)~=0)
            disp('Parameters too close to each other, drawing new ones.');
        end
        Wt = rand(Npar,Nmodes)*4-2;
    end
    modeT = mod(randperm(N),Nmodes)+1;
    for i=1:N
        Y(i) = X(i,:)*Wt(:,modeT(i));
    end

    Noise = randn(N,1)*sigmaNoise;

    Y = Y+Noise;
    
    % klinreg
    
    tic; 

    [Wtmp lambda]=klinreg(X,Y,Nmodes);

    lambda_ordered=zeros(N,1);
    for p=1:Nmodes
        
        [mm j] = min(diag((Wtmp-Wt(:,p)*ones(1,Nmodes-p+1))'*(Wtmp-Wt(:,p)*ones(1,Nmodes-p+1))));
        Wtmp_ordered(:,p) = Wtmp(:,j);        
        Wtmp(:,j)=[];
        imp = find(lambda==j);
        lambda_ordered(imp) = p;
        impp=find(lambda>j);
        lambda(imp)=0;
        lambda(impp)=lambda(impp)-1;
    end
    % Computing time
    comp_time(expe)=toc; 
    % Normalized parametric error 
    NMSEparam(expe) = sum(diag((Wtmp_ordered-Wt)'*(Wtmp_ordered-Wt))./diag(Wt'*Wt))/Nmodes;
    % Classification error
    ClassifErrRate(expe) = sum(lambda_ordered~=modeT')/N;
	   
	if DO_OTHERS
	    % ME-MCS
	    tic;
	    
	     [Wtmp lambda]=me_mcs(X,Y,Nmodes);
	    lambda_ordered=zeros(N,1);
	    for p=1:Nmodes
		
		[mm j] = min(diag((Wtmp-Wt(:,p)*ones(1,Nmodes-p+1))'*(Wtmp-Wt(:,p)*ones(1,Nmodes-p+1))));
		Wtmp_ordered(:,p) = Wtmp(:,j);        
		Wtmp(:,j)=[];
		imp = find(lambda==j);
		lambda_ordered(imp) = p;
		impp=find(lambda>j);
		lambda(imp)=0;
		lambda(impp)=lambda(impp)-1;
	    end
	    comp_timeME(expe)=toc;
	    NMSEparamME(expe) = sum(diag((Wtmp_ordered-Wt)'*(Wtmp_ordered-Wt))./diag(Wt'*Wt))/Nmodes;
	    ClassifErrRateME(expe) = sum(lambda_ordered~=modeT')/N;
	   
	   % PE-MCS
	    tic;
	    
	     [Wtmp lambda]=pe_mcs(X,Y,Nmodes);
	    lambda_ordered=zeros(N,1);
	    for p=1:Nmodes
		
		[mm j] = min(diag((Wtmp-Wt(:,p)*ones(1,Nmodes-p+1))'*(Wtmp-Wt(:,p)*ones(1,Nmodes-p+1))));
		Wtmp_ordered(:,p) = Wtmp(:,j);        
		Wtmp(:,j)=[];
		imp = find(lambda==j);
		lambda_ordered(imp) = p;
		impp=find(lambda>j);
		lambda(imp)=0;
		lambda(impp)=lambda(impp)-1;
	    end
	    comp_timePE(expe)=toc;
	    NMSEparamPE(expe) = sum(diag((Wtmp_ordered-Wt)'*(Wtmp_ordered-Wt))./diag(Wt'*Wt))/Nmodes;
	    ClassifErrRatePE(expe) = sum(lambda_ordered~=modeT')/N;
	end
    
    %% Reference model
    lambda_ordered=modeT;
    tic; 
    for p=1:Nmodes
        Wtmp_ordered(:,p) = pinv(X(modeT==p,:))*Y(modeT==p);        
    end
    [eRef lambda_ordered]=min((Y*ones(1,Nmodes)-X*Wtmp_ordered).^2, [],2);
    comp_timeRef(expe)=toc;
    NMSEparamRef(expe) = sum(diag((Wtmp_ordered-Wt)'*(Wtmp_ordered-Wt))./diag(Wt'*Wt))/Nmodes;
    ClassifErrRateRef(expe) = sum(lambda_ordered~=modeT')/N;
    
end

% Results of the reference model
MeanTimeRef = mean(comp_timeRef)
StdTimeRef = std(comp_timeRef)
MaxTimeRef = max(comp_timeRef)
MeanNMSERef = mean(NMSEparamRef)
StdNMSERef = std(NMSEparam-MeanNMSERef)
MaxNMSERef = max(NMSEparamRef)
MeanClassifErrRatePourcentRef = 100*mean(ClassifErrRateRef);
StdClassifErrRatePourcentRef = std(100*ClassifErrRateRef-MeanClassifErrRatePourcentRef);
MaxClassifErrRatePourcentRef = 100*max(ClassifErrRateRef);

% Results of k-linreg
MeanTime = mean(comp_time)
StdTime = std(comp_time)
MaxTime = max(comp_time)
MeanNMSE = mean(NMSEparam)
StdNMSE = std(NMSEparam-MeanNMSE)
MaxNMSE = max(NMSEparam)
MeanClassifErrRatePourcent = 100*mean(ClassifErrRate);
StdClassifErrRatePourcent = std(100*ClassifErrRate-MeanClassifErrRatePourcent);
MaxClassifErrRatePourcent = 100*max(ClassifErrRate);

if DO_OTHERS
	% Results of ME-MCS
	MeanTimeME = mean(comp_timeME)
	StdTimeME = std(comp_timeME)
	MaxTimeME = max(comp_timeME)
	MeanNMSEME = mean(NMSEparamME)
	StdNMSEME = std(NMSEparamME-MeanNMSEME)
	MaxNMSEME = max(NMSEparamME)
	MeanClassifErrRatePourcentME = 100*mean(ClassifErrRateME);
	StdClassifErrRatePourcentME = std(100*ClassifErrRateME-MeanClassifErrRatePourcentME);
	MaxClassifErrRatePourcentME = 100*max(ClassifErrRateME);

	% Results of PE-MCS
	MeanTimePE = mean(comp_timePE)
	StdTimePE = std(comp_timePE)
	MaxTimePE = max(comp_timePE)
	MeanNMSEPE = mean(NMSEparamPE)
	StdNMSEPE = std(NMSEparamPE-MeanNMSEPE)
	MaxNMSEPE = max(NMSEparamPE)
	MeanClassifErrRatePourcentPE = 100*mean(ClassifErrRatePE);
	StdClassifErrRatePourcentPE = std(100*ClassifErrRatePE-MeanClassifErrRatePourcentPE);
	MaxClassifErrRatePourcentPE = 100*max(ClassifErrRatePE);
end
