% Experiments of section 5.2: High dimensions
%
% Measure the computing time of klinreg over 100 randomly generated problems for different dimensions
%
% You must have klinreg_matlab installed and in the PATH
% to run this test.
%
clear all

N = 10000;

% List of values for the dimension p
NparT = [3,5,10,20,30,50,75,100,150,200]; 

% List of values for the number of mode n
NmodesT = [2,3,5];

Nexpe=100; % number of experiments for each pair (n,p)

% Prepare storage for results
% BigXXX is a matrix of the values of the quantity XXX
% for all values of n (row-wise) and p (column-wise)
BigTime = zeros(length(NmodesT),length(NparT));
BigTimeStd = zeros(length(NmodesT),length(NparT));
BigNMSE = zeros(length(NmodesT),length(NparT));
BigNMSEStd = zeros(length(NmodesT),length(NparT));

% Parameters of klinreg_simple
restarts = zeros(length(NmodesT),length(NparT));
wmax = 100;

Nmi = 0;
for Nmodes=NmodesT	% for each number of modes
    Nmi = Nmi+1
    Npari = 0;
    for Npar=NparT	% for each dimension
        Npari = Npari + 1;

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

	ClassifErrRate = zeros(Nexpe,1);
	Time = zeros(Nexpe,1);

	% Compute optimal number of restarts
	restarts(Nmi,Npari) = rstar(0.001, N,Nmodes,Npar);

	% Do Nexpe experiments to compute averages
	for expe=1:Nexpe
	    expe
	    
	    % data
	    X = rand(N,Npar)*2-1; 
	    
	    % true parameters
	    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

	    % mode
	    modeT = mod(randperm(N),Nmodes)+1;
	    for i=1:N
		Y(i) = X(i,:)*Wt(:,modeT(i));
	    end
	    varianceNoise = 0.01; 
	    Noise = randn(N,1)*sqrt(varianceNoise);
	    Y = Y+Noise;

	    % k-linreg
	    
	    tic; 
	    [Wtmp, lambda] = klinreg_simple(X,Y,Nmodes,restarts(Nmi,Npari),wmax);
	    
	    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
	    Time(expe)=toc;
	    NMSE(expe) = sum(diag((Wtmp_ordered-Wt)'*(Wtmp_ordered-Wt))./diag(Wt'*Wt))/Nmodes;
	    ClassifErrRate(expe) = sum(lambda_ordered~=modeT')/N;
	end

	NMSEthresh = NMSE;
	NMSEthresh(NMSE>1)=[];
	Nfails(Nmi,Npari) = length(NMSE)-length(NMSEthresh);
	NMSE = NMSEthresh;


	MeanTime = mean(Time);
	StdTime = std(Time);
	MaxTime = max(Time);
	MeanNMSE = mean(NMSE);
	StdNMSE = std(NMSE-MeanNMSE);
	MaxNMSE = max(NMSE);

	% Store the results 
	BigTime(Nmi,Npari) = MeanTime; 	
	BigTimeStd(Nmi,Npari) = StdTime;
	BigNMSE(Nmi,Npari) = MeanNMSE;
	BigNMSEStd(Nmi,Npari) = StdNMSE;    
    end
end
disp('Computing time of k-LinReg for different values of n (column-wise) and p (row-wise)');
CpuTime_vs_n_p = BigTime
n_for_each_row_is = NmodesT
p_for_each_column_is = NparT
