% It is based on the code of the file compareProbVsClassicBernoulli.m, with Copyright Sohan Seth sohan.seth@hiit.fi 

% PUT bin R on on the Matlab path


close all,  
saveFigure = 0; plotFigure = 1; maxTrial = 10; % 10 for figure, 100 for table with plotFigure = 0;
maxStart = 10; % Run PAA maxStart times and take best solution

computeCostBernoulli = @(X, P)(sum(sum(-X.*log(double(P + eps)))) ...
    + sum(sum(-(1 - X).*log(1 - double(P) + eps))) );
computeCostPoisson = @(X, R)(sum(sum(- X .* log(R + eps) + R)));
computeCostStochastic = @(X, P)...
    (- sum(sum(X .* log(P + eps),2)));

%% Bernoulli data, multiple data multiple trials
%setenv('DYLD_LIBRARY_PATH', '/usr/local/bin/'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%rng default, %version higher matlab
%rand('twister',0) %matlab version lower
rng(1, 'twister');
%https://www.mathworks.com/help/matlab/math/generate-random-numbers-that-are-repeatable.html

clear nFeatSam
d = 10;         % Dimensionality
K = 6;          % Number of archetypes
n = 100;        % Number of observations
gammaParam = 0.8; % Concentration of Dirichlet distribution for H
options = generate_options();
options.verbose = false; options.maxIter = 10000;
FONTSIZE = 6;
if plotFigure; hMain = figure('papersize',[8.5 11],'paperposition',[0 0 8 2]); %colormap(flipud(gray)); 
colormap gray
end

errpaa=[];
erraanos=[];
erraans=[];
for countTrial = 1:maxTrial
    matFeatLat = rand(d, K) > 0.7; %archetypes
     matFeatLatn = imnoise(double(matFeatLat),'salt & pepper', 0.05); %salt an peper noise
    %alphas
    matLatSam = gamrnd(gammaParam, 1, K, n); matLatSam = bsxfun(@rdivide, matLatSam, sum(matLatSam));
   
    matFeatSam = imnoise(matFeatLatn * matLatSam,'gaussian') >0.5;
    
    if plotFigure
        subplot(4,maxTrial,countTrial)
        imagesc(matFeatLat)
        if countTrial == 1
            ylabel('True','fontsize',FONTSIZE)
        end
        title(num2str(countTrial),'fontsize',FONTSIZE)
        set(gca,'fontsize',FONTSIZE,'xtick',[],'ytick',[])
    end
    
    % Run PAA and take the best solution
    matSamLatEM = cell(maxStart,1);
    matLatSamEM = cell(maxStart,1);
    objEM = cell(maxStart,1);
    for countEM = 1:maxStart
        [matSamLatEM{countEM}, matLatSamEM{countEM}, objEM{countEM}] = paa_Bernoulli(matFeatSam, K, options);
    end
    [~, ind] = min(cellfun(@(x)(x(end)), objEM));
    matSamLatEM = matSamLatEM{ind};
    matLatSamEM = matLatSamEM{ind};
    
   

	[TTT,err]=permmatch(double(matFeatSam * matSamLatEM),double(matFeatLat),K);  
 
    errpaa=[errpaa err];
    
    if plotFigure
        subplot(4,maxTrial,maxTrial + countTrial)
        imagesc(TTT)% > 0.5)
      %  text(find(sum((bsxfun(@eq,ind,ind'))) ~= 1),0*double(find(sum((bsxfun(@eq,ind,ind'))) ~= 1) > 0),'o','fontsize',FONTSIZE/2)
        if countTrial == 1
            ylabel('PAA','fontsize',FONTSIZE)
        end
        set(gca,'fontsize',FONTSIZE,'xtick',[],'ytick',[])
    end
    
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%AA without standarizing data, with RAW data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[matSamLatClassicnos, matLatSamClassicnos, objnos] = nonstand_aa(matFeatSam, K);

 [TTT,err]=permmatch(double(matFeatSam * matSamLatClassicnos),double(matFeatLat),K);  
 
    erraanos=[erraanos err];
   

    
    if plotFigure
        subplot(4,maxTrial,2*maxTrial + countTrial)
        imagesc(TTT)% > 0.5)
        text(find(sum((bsxfun(@eq,ind,ind'))) ~= 1),0*double(find(sum((bsxfun(@eq,ind,ind'))) ~= 1) > 0),'o','fontsize',FONTSIZE/2)
        if countTrial == 1
            ylabel('AA','fontsize',FONTSIZE)
        end
        set(gca,'fontsize',FONTSIZE,'xtick',[],'ytick',[])
    end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Archetypoids RAW data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[matSamLataans, matLatSamaans, objaans] = ada3(matFeatSam, K);




 [TTT,err]=permmatch(double(matFeatSam(:,matSamLataans)),double(matFeatLat),K);  
 
    erraans=[erraans err];
   
  
    
    if plotFigure
        subplot(4,maxTrial,3*maxTrial + countTrial)
        imagesc(TTT)% > 0.5)
        text(find(sum((bsxfun(@eq,ind,ind'))) ~= 1),0*double(find(sum((bsxfun(@eq,ind,ind'))) ~= 1) > 0),'o','fontsize',FONTSIZE/2)
        if countTrial == 1
            ylabel('ADA ','fontsize',FONTSIZE)
        end
        set(gca,'fontsize',FONTSIZE,'xtick',[],'ytick',[])
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


e=[ errpaa; erraanos; erraans]

mean(e')
sqrt(var(e'))

median(e')

%10 first

%e =

%     3     2     3     4     5     5     4     4     3     2
%     3     2     2     4     4     2     3     4     2     2
%     3     2     2     3     2     2     3     4     2     2


%ans =

 %   3.5000    2.8000    2.5000


%ans =

 %   1.0801    0.9189    0.7071


%ans =

 %   3.5000    2.5000    2.0000

%100 first


%e =

  


%  Columns 1 through 21

 %    3     2     3     4     5     5     4     4     3     2     4     3     5     5     6     5     5     2     5     5     2
 %    3     2     2     4     4     2     3     4     2     2     4     3     5     5     5     5     5     3     2     5     4
 %    3     2     2     3     2     2     3     4     2     2     3     3     1     5     2     3     5     2     3     4     1

 % Columns 22 through 42

  %   7     4     4     6     5     2     2     6     4     4     4     2     7     5     6     2     7     1     4     2     4
  %   8     1     7     4     5     2     3     3     7     4     5     2     6     5     6     1     7     0     4     2     4
  %   4     1     3     3     2     2     3     3     5     3     3     3     6     5     6     0     7     0     4     2     3

 % Columns 43 through 63

  %   7     5     4     6     2     4     2     4     2     2     4     2     6     5     2     3     7     2     3     5     5
  %   2     4     4     6     2     2     2     3     2     4     3     2     3     2     3     1     5     1     6     4     5
  %   2     4     4     7     2     2     2     2     2     2     3     2     6     2     3     1     5     2     3     3     7

  %Columns 64 through 84

   %  7     5     3     3    10     3     5     6     2     3     5     7     2     6     2     3     6     2     3     5     4
    % 8     3     1     2     8     0     2     5     1     3     5     7     2     6     2     3     5     2     2     5     4
    % 8     5     2     2     7     0     2     4     2     4     4     5     2     6     1     3     4     1     2     5     4

  %Columns 85 through 100

   %  8     5     4     2     2     2     3    10     4     8     5     6     3     6     3     5
    % 7     4     2     2     1     1     2    10     2     5     1     6     3     6     2     3
    % 7     5     1     2     1     1     2    10     2     5     1     6     2     5     2     3

%ans =

 %   4.2000    3.5900    3.1900


%ans =

%    1.8641    1.9904    1.8786


%ans =

 %    4     3     3

