%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example script for the demonstartion of Functional Embedding (FEM)
% Data: Tecator data (Borggaard and Thodberg, 1992), available at
%       http://lib.stat.cmu.edu/datasets/tecator
%       For demo. purpose, please DL "spectrometric.dat" from the page 
%       where you find this script
% Required soft/package:
%       1. MATLAB 
%       2. PACE: http://anson.ucdavis.edu/~mueller/data/programs.html
% Submission mode: after DL data and PACE, 
%                  remodify MATLAB path list 
%                  (for know-how see PACE website)
%       1. Interactive: (a)loading the data 
%                       (b)highlight block (steps) of codes as you
%                          prefer, press F9 (under window OS)
%       2. Batch: F5 (under window OS)                                         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear, close all;
% (a)
tecator=load('spectrometric.dat');
x=tecator(:,1:100); 
% Spectrometric data: intensity of 215 subj.(row) over 100 chanels(col).
% note that the first 172 rows are labeled as training set.
%           the last  43  rows are labeled as test     set.
y=tecator(:,101);   % Response: fat content

% (b-1) Creation of groups: gp=1 if y_i> median(y); 0, otherwise. 
medy=median(y);
gp=zeros(size(y));
gp(y>=medy)=1;

% (b-2) Random scrambling of the ordering of 100 chanels.
rand('state',1001);
[tmp,perindx]=sort(rand([1,100]));
newx=x(:,perindx); % Scrambled data

% (b-3):FEM step
%  utilizing classical MDS to extract the univariate 
%  configuration on dissimilarity matrix: I-R 
Tot_d=pdist(newx','correlation');
[u,tmp]=cmdscale(Tot_d);
u=u(:,1);
[tempu,sindx]=sort(u);
% For better display, reverse the orientation. (optional)
sindx=sindx(100:-1:1);

% Reshuffle the 
snewx=newx(:,sindx); % Scrambled data after cleaning by FEM
su=linspace(1,100,100); % Indices for the chanels.

figure(101)
subplot(3,1,1)
plot(1:100,x','k-');
axis([0,100,1.5,6.2])
title('Original data')
subplot(3,1,2)
plot(1:100,newx','k-');
axis([0,100,1.5,6.2])
title('Scrambled data')
subplot(3,1,3)
plot(1:100,snewx','k-');
axis([0,100,1.5,6.2])
title('Data cleaned after FEM')

% Reference: 
% 1. Wu, P.S., and Mueller, H.G. (2009) "Functional Embedding]{Functional Embedding for the Classification of
%    Gene Expression Profiles" submitted to Bioinformatics
% 2. PACE. 
