%% DON"T USE is ROB's MAGIC STUFF
%%restoredefaultpath; matlabrc 
clear all
clc
close all


disp('Startup from LOeKS directory');
cd 'E:\MATL_DIR\ANALYSIS_LOEK_HEIJNEN\LINVECMOD\'
addpath( ...
    genpath(fullfile('E:\','MATL_DIR\ANALYSIS_LOEK_HEIJNEN\LINVECMOD\')));
cd  'E:\MATL_DIR\ANALYSIS_LOEK_HEIJNEN\LINVECMOD\'

%------- CALCULATE ALL PERMUTATIONS -------%
clc
clear all
close all


%sig = 1.5; %%% .75 is a 9 unit gabor patch
%sizes       = 18;
% pos         = 6:1:(sizes-(6-1));
% MIN_DIST    = 5;
% height      = [0.0 linspace(0.1,1.0,10)]; %%% alpha

sig = 2.5; 
sizes = 20;


pos         = 6:1:(sizes-(6-1));
MIN_DIST    = 5;
height      = [0.0 linspace(0.1,1.0,6)]; %%% alpha


T1=[];
T2=[];


%%% Caculate all position combinations of T1 and T2
for i=1:1:length(pos)
            T1 = [ T1   repmat(pos(i), [1 , length(pos)]) ];       
            T2 = [ T2  pos];
end
T1T2= [T1' T2'];

%%% Make sure the top of T1 does not ovelap that of T2
deltaT=(T2'-T1');
dist=abs(deltaT);
T1T2=T1T2(dist >= MIN_DIST,:);

%%% caculate all combinations of T1 and T2 when combined with all alphas
N=length(T1T2);
B = height(ones(1,N),:);
C = reshape(B, numel(B), 1); % gives a1 a1 a2 a2 a3 a3
%C = reshape(B', numel(B), 1); % gives a1 a2 a3 a1 a2 a3
alpha = C;

T1T2= repmat(T1T2,length(height),1);
T1T2=[T1T2 alpha];

T1p      =   T1T2(:,1); 
T2p      =   T1T2(:,2); 
deltaT   =   (T2p - T1p);
alpha    =   T1T2(:,3); 
avgTp    =   T1p +  alpha.* deltaT;
difTp    =   (1 - alpha).* deltaT;


figure(44);
cla

subplot(231)
plot(T1p(:,1),'.');
title([ 'T1 INPUT']);

subplot(232)
plot(T2p(:,1),'.');
title([ 'T2 INPUT']);

subplot(233)
plot(deltaT(:,1),'.');
title([ 'deltaT']);

subplot(235)
plot(avgTp,'.');
title([ 'avgT']);

subplot(236)
plot(difTp,'.');
title([ 'difT']);

%%
%------- GENERATE INPUT / OUTPUT DATA -------%

T1m=zeros(sizes);
T2m=zeros(sizes);
Tavgm=zeros(sizes);
Tdifm=zeros(sizes);

samples_length = length(T1T2);
p = zeros(sizes,samples_length);
p1 = zeros(sizes,samples_length);
p2 = zeros(sizes,samples_length);
tavg = zeros(sizes,samples_length);
tdif = zeros(sizes,samples_length);

for k=1:1:length(T1T2)
        T1m =                ONE_D_gabor( sizes, sig, T1p(k) );
        T2m = alpha(k)   *   ONE_D_gabor( sizes, sig, T2p(k) );

        Tavgm   =           ONE_D_gabor( sizes, sig, avgTp(k) );
        
        x_shift =           abs(min(difTp))+min(pos)-4;%-2; %%%% make sure Tdifm postion is not negative
        Tdifm   =           ONE_D_gabor( sizes, sig, difTp(k) + x_shift);

          % INPUT %
          p(:,k)  = reshape(T1m + T2m, numel(T1m + T2m), 1);  
          p1(:,k) = reshape(T1m, numel(T1m), 1);  
          p2(:,k) = reshape(T2m, numel(T2m), 1);  

          % OUTPUT %
          tavg(:,k)= reshape(Tavgm, numel(Tavgm), 1); 
          tdif(:,k)= reshape(Tdifm, numel(Tdifm), 1);  
end

p=[p1; p2];
t=[ tavg; tdif];

[pn,stdp] = mapminmax(p);
[tn,stdt] = mapminmax(t);

    net = newff(minmax(pn),[length(t(:,1))],{'tansig'},'trainbr');

     net.trainParam.show = 3;
     net.trainParam.epochs = 20;
     net.trainParam.goal = 10/10^20;
     net.trainParam.mem_reduc =1;
     %net.trainParam.lr= 0.0001;
     %%net.divideFcn = '';
    net.performFcn = 'msereg';
    net.performParam.ratio = 2;

% net.trainParam.mu= 1.5 %Initial Mu
% net.trainParam.mu_dec = .9 %Mu decrease factor
% net.trainParam.mu_inc = 2 %Mu increase factor
% net.trainParam.mu_max = 1e10 %Maximum Mu

    
    %-- Train the network --%
    randn('seed',192736547);
    net = init(net);
                                           %VAL   %TEST
    [TrainV,ValV,TestV] = dividevec(pn,tn,25/100,25/100);
    
    [net,tr] = train(net,TrainV.P,TrainV.T,[],[],ValV,TestV);

    %-- Simulate the network --%
    [A1,Pf,Af,E,info.train.performance]         = sim(net,TrainV.P, [],[],TrainV.T);
    [A2,Pf,Af,E,info.validation.performance]    = sim(net,ValV.P,   [],[],ValV.T);
    [A3,Pf,Af,E,info.test.performance]          = sim(net,TestV.P,  [],[],TestV.T);

    Y1 = mapminmax('reverse',A1,stdt);
    Y2 = mapminmax('reverse',A2,stdt);
    Y3 = mapminmax('reverse',A3,stdt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(15)
cla
    tP=t(:,TrainV.indices); %% teacher data
    oP=Y1;                  %% output simulation
    Q=numel(oP)/2;
            subplot(321)
            postregV2(oP(1:Q),tP(1:Q));
            title('Tavg')
            axis square
            grid
            subplot(322)
            postregV1(oP(Q+1:Q*2),tP(Q+1:Q*2));
            title('Tdif')
            axis square
            grid
        
    tV=t(:,ValV.indices);   %% teacher data
    oV=Y2;                  %% output simulation
    Q=numel(oV)/2;
          subplot(323)
            postregV2(oV(1:Q),tV(1:Q));
            title('Tavg')
            axis square
            grid
          subplot(324)
            postregV2(oV(Q+1:Q*2),tV(Q+1:Q*2));
            title('Tdif')
            axis square
            grid
   
    tT=t(:,TestV.indices);  %% teacher data
    oT=Y3;                  %% output simulation
    Q=numel(oT)/2;
          subplot(325)
            postregV1(oT(1:Q),tT(1:Q));
            title('Tavg')
            axis square
            grid

          subplot(326)
            postregV1(oT(Q+1:Q*2),tT(Q+1:Q*2));
            title('Tdif')
            axis square
            grid     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%%
    %-- SIMULATION TRAINED DATA--%
clc;
figure(18);
hold off;

for k=1:1:length(oP)
cla  

pt1 = p1(:,TrainV.indices);
pt2 = p2(:,TrainV.indices);
q = TrainV.indices(k);

subplot(311)

    plot(pt1(:,k),'k-');
    hold on;
    plot(pt2(:,k),'b-');
    xlim([0.5 sizes+0.5]);
    ylim([-0.3 1.3]);
    title([ 'INPUT T1: ' num2str(T1p(q))  ' & T2: ' num2str(T2p(q))  ' ALPHA: ' num2str(alpha(q))]);
    hold off;
    
subplot(323)
    plot(oP(1:sizes,k),'r-');
    hold on;
    plot(tP(1:sizes,k),'g-');
    xlim([0.5 sizes+0.5]);
    ylim([-0.3 1.3]);
    title([ 'Tavg: '  num2str(avgTp(q)) ]);
    hold off;
    
subplot(324)
    plot(oP(sizes + 1 : sizes*2,k),'r-');
    hold on;
    plot(tP(sizes+1:sizes*2,k),'g-');
    xlim([0.5 sizes+0.5]);
    ylim([-0.3 1.3]);
    title([ 'Tdif: '  num2str( difTp(q)  ) ]);%%% + abs(min(difTp)) + min(pos)-2
 
    Q=1:2:18;
    set(gca,'XTick',Q);
    set(gca,'xTickLabel',Q - x_shift);
    hold off;
    
 subplot(325)
    postregV1(oP(1:sizes,k)',tP(1:sizes,k)')

    
 subplot(326)
    postregV1(oP(sizes+1:sizes*2,k)',tP(sizes+1:sizes*2,k)')

pause(0.3)

end

