
close all

%%%%%%%%%%%%%%%%%%%%
% Define variables %
%%%%%%%%%%%%%%%%%%%%

selector = 1; % 1 for sources, 2 for mixtures


rows = 2000;
columns = 100;
sets = 2;
sources = rand(rows, columns);
mixtures = rand(rows,columns);

dat = zeros(rows, columns, sets);

figure
imagesc(sources)


%%%%%%%%%%%
% Sources %
%%%%%%%%%%%

% Increase weighting of sine feature and check cFLs


for i = 1:rows
   x = [0.1:0.1:10];
   n = rand(1,100);
   y = x + n;
   
   
   if i < 500
        z = sin(y);
          
   elseif i >= 500 && i < 1000
       z = 4*sin(y);
       
   elseif i >= 1000 && i < 1500
       z = 8*sin(y);
       
   elseif i >= 1500 && i <= 2000
       z = 16*sin(y);
          
   end
   
   sources(i,:) = sources(i,:) + z;
end
   
%%%%%%%%%%%%
% Mixtures %
%%%%%%%%%%%%

temps = zeros(rows,1);

for i = 1:rows
   x = [0.1:0.1:10];
   n = rand(1,100);
   y = x + n;
   
   T = 10*(1/50);
   fs = 500;
   t = 0:1/fs:T-1/fs;
   st = sawtooth(2*pi*50*t);
   st = st + n;
   temp = round(rand(1)) + 1;
   


   if i < 500
       
       z = sin(y) + 6*temp*st; %+ step;
          
   elseif i >= 500 && i < 1000
       z = 4*sin(y) + 6*temp*st; %+ step;
       
   elseif i >= 1000 && i < 1500
       z = 8*sin(y) + 6*temp*st; %+ step;
       
   elseif i >= 1500 && i <= 2000
       z = 16*sin(y) + 6*temp*st; %+ step;
          
   end
   
   mixtures(i,:) = mixtures(i,:) + z;
   temps(i,1) = round(temp);
end
   


%%%%%%%%
% Data %
%%%%%%%%

dat(:,:,1) = sources;
dat(:,:,2) = mixtures;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following uses the SVDS function to compute solutions for specific  %
% values of k. This decomposition is applied to the whole data matrix and %
% yields k eigenvalues (in Va), k un-normalized eigenvectors (V2), and    %
% component loadings for each row of the original data matrix (cTCs).     %         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
% plot eigenvalues %
%%%%%%%%%%%%%%%%%%%%

dat = dat(:,:,selector);
data = zscore(dat);

k1 = size(data,1);
opts = struct('disp', 0, 'maxit', 1000);

[Ua1, La1, Va1] = svds(data, k1, 'L', opts);

La1 = La1 ./ sqrt(size(data, 1) - 1); % corrects for differences in the rank of the data matrix
La1 = La1.^2; % See https://en.wikipedia.org/wiki/Singular_value_decomposition; diagonal is the square root of the eigenvalues
V2a = Va1*sqrt(La1); % Un-normalizes the eigenvector matrix
cFLs = data(:,:)*V2a(:,1:4); 
LGD1 = diag(La1); % Vector of eigenvalues for scree plot, AOC estimation

eigenvectors = V2a;
eigenvalues = LGD1;

%%%%%%%%%
% Plots %
%%%%%%%%%


males = find(temps==1);
females = find(temps==2);

figure
plot(squeeze(dat(males(1),:)))
hold on
plot(squeeze(dat(females(1),:)))

%figure
plot(squeeze(dat(males(end),:)))
%hold on
plot(squeeze(dat(females(end),:)))


figure
imagesc(squeeze(dat))


figure
plot(eigenvalues(1:5))

figure
plot(eigenvectors(:,1))
hold on
plot(eigenvectors(:,2))

figure
plot(cFLs)

figure
plot(st)
hold on
plot(sin(y))

%%

%---------------%
% Design Matrix %
%---------------%

DM = zeros(2000,4); % ID var, Sex, age, SES


for i = 1:2000
    
    if i < 500
        DM(i,1) = 2 + rand(1);
        DM(i,2) = 1 + temps(i,1);
        DM(i,3) = 100*rand(1);
        DM(i,4) = 3*rand(1);
        
    elseif i >= 500 && i < 1000
        DM(i,1) = 3 + rand(1);
        DM(i,2) = 1 + temps(i,1);
        DM(i,3) = 100*rand(1);
        DM(i,4) = 3*rand(1);
        
    elseif i >= 1000 && i < 1500
        DM(i,1) = 4 + rand(1);
        DM(i,2) = 1 + temps(i,1);
        DM(i,3) = 100*rand(1);
        DM(i,4) = 3*rand(1);
        
    elseif i >=1500 && i <=2000
        DM(i,1) = 5 + rand(1);
        DM(i,2) = 1 + temps(i,1);
        DM(i,3) = 100*rand(1);
        DM(i,4) = 3*rand(1);
        
    end
end
  
figure
imagesc(DM, [0,10])


lm_output = cell(1,4);
for i = 1:size(DM,2)
    
    FL = cFLs(:,i);
    Disease = DM(:,1);
    Sex = DM(:,2);
    Age = DM(:,3);
    SES = DM(:,4);
    tbl = table(Disease,Sex,Age,SES,FL);  
    % Fit Linear Model
    lm = fitlm(tbl,'FL~Disease+Sex+Age+SES');
    lm_output{1,i} = lm;
end
    