%% Example
%% Load diabetes data
clear;

    data='EEG';
    goal='Autoenc';
    k=10;
    ridgeParams=10;
    shiftEMG=0.05;
    shiftEEG=0.2;
    eegOffset=1;

    load '/home/jph/Uni/masterarbeit/messor/matlabData/MB2200msWindowEMG1sWindowEEG50msShiftEMG200msShiftEEGFreq0to200Pause0pBurg50.mat';

    clear classification;
    eeg=false;
    if strcmp(data,'EEG')
        eeg=true;
        trainingData=trainingDataEEG(eegOffset+1:end,:);
    elseif strcmp(data,'EMG')
        trainingData=trainingDataEMG;
    elseif strcmp(data,'LF')
        eeg=true;
        trainingData=trainingDataEEGlf(eegOffset+1:end,:);
    else
        error('only EEG, EMG and LF are valid inputs for data');
    end
    
    factor=size(trainingDataEMG,1)/(size(trainingData,1)+eegOffset);

    clear trainingDataEEG;
    clear trainingDataEEGlf;
    clear trainingDataEMG;

    if strcmp(goal,'Autoenc')
        predicted=shiftingMean(synergiesAutoenc,factor);
    elseif strcmp(goal,'PCA')
        predicted=shiftingMean(synergiesPCA,factor);
    elseif strcmp(goal,'NNMF')
        predicted=shiftingMean(synergiesNNMF,factor);
    elseif strcmp(goal,'kin')
        predicted=shiftingMean(kinematics,factor);
    else
        error('only kin, Autoenc, PCA and NNMF are valid inputs for goal');
    end

    clear kinematics;
    clear synergies*;

    if eeg
        predicted=predicted(1:end-eegOffset,:);
    end

X = trainingData;
y = predicted;

n = size(X,1);

%% Augment the design matrix with large number of additional variables
Xn = [X, randn(n, 4e3)];
p  = size(Xn,2);

fprintf('*****************************************************************************\n');
fprintf('* fastridge example using diabetes data (augmented with extra noise variables):\n');
fprintf('  - n = %d\n', n);
fprintf('  - p = %d\n', p);
fprintf('\n');

%% Produce a regularisation path
fprintf('* Producing a regularisation path of 100 betas using fastridge ');
tic;
[beta, b0, tau2, DOF, lambda, score] = fastridge(Xn, y)%, 'path', 1e2);
fprintf('(%.2fs)\n', toc);

[~,I] = min(score);
fprintf('  - DOF of best model    = %.2f\n', DOF(I));
fprintf('  - lambda of best model = %.2f\n', lambda(I));

displaypath(beta, DOF, score);

fprintf('\n');

%% Fit model with specified DOF
fprintf('* Fitting model with %.2f DOF ', DOF(I));
tic;
[b, ~, ~, ~, ~, ~] = fastridge(Xn, y, 'DOF', DOF(I));
fprintf('(%.2fs)\n', toc);
fprintf('  - L2 norm between this model and best model from path = %.2f\n', norm(b-beta(:,I)));

fprintf('\n');

%% Fit model with specified lambda
fprintf('* Fitting model with lambda=%.2f ', lambda(I));
tic;
[b, ~, ~, ~, ~, ~] = fastridge(Xn, y, 'lambda', lambda(I));
fprintf('(%.2fs)\n', toc);
fprintf('  - L2 norm between this model and best model from path = %.2f\n', norm(b-beta(:,I)));

fprintf('\n');

%% Search for best model using MML criterion
fprintf('* Fitting model using MML (searching for lambda) ');
tic;
[b_srch, ~, ~, DOF_srch, lambda_srch] = fastridge(Xn, y);
fprintf('(%.2fs)\n', toc);
fprintf('  - DOF of this model    = %.2f\n', DOF_srch);
fprintf('  - lambda of this model = %.2f\n', lambda_srch);
fprintf('  - L2 norm between this model and best model from path = %.2f\n', norm(b_srch-beta(:,I)));

fprintf('\n');

%% Search for best model using AICc criterion
fprintf('* Fitting model using AICc (searching for lambda) ');
tic;
[b_srch, ~, ~, DOF_srch, lambda_srch] = fastridge(Xn, y, 'criterion', 'aicc');
fprintf('(%.2fs)\n', toc);
fprintf('  - DOF of this model    = %.2f\n', DOF_srch);
fprintf('  - lambda of this model = %.2f\n', lambda_srch);
fprintf('  - L2 norm between this model and best model from MML path = %.2f\n', norm(b_srch-beta(:,I)));

fprintf('\n');

%% Finally fit a model using the standard ridge
fprintf('* Fitting model with lambda=%.2f using MATLAB Ridge ', lambda(I));
tic;
b_matlab = ridge(y,Xn,lambda(I),0);
fprintf('(%.2fs)\n', toc);
fprintf('  - L2 norm between this model and best model from path = %.2f\n', norm(b_matlab(2:end)-beta(:,I)));