%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Experiment with gradient ascent                                         %
% Project: sub-euclidean distance for music similarity                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load features
%feature_file = 'rel_music_raw_features.mat';
feature_file = 'rel_music_raw_features+simdata_ISMIR12.mat';

vars = whos('-file', feature_file);
A = load(feature_file,vars(1).name,vars(2).name,vars(3).name,vars(4).name);
raw_features = A.(vars(1).name);
indices      = A.(vars(2).name);
tst_inx      = A.(vars(3).name);
trn_inx      = A.(vars(4).name);
%% Params setting
dmr    = [0 5 10 20 30 50];    % dimension reduction by PCA
ws     = [0 5 10 20 30 50 70]; % window size
% parameters of rbm (if it is used for extraction)
hidNum = 0;
lr_1   = 0;
lr_2   = 0;
mmt    = 0;
cost   = 0;
%% Select parameters (if grid-search is not applied)
 di  = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If grid search is define
% Define directory to save parameters & results
if ~isempty(findstr('WIN',computer()))
        dir    = 'C:\Pros\Experiments\ISMIR_2013\grad\'; % In windows platform
        dlm    = '\';
elseif ~isempty(findstr('linux',computer())) || ~isempty(findstr('LNX',computer()))
        dir    = '/home/funzi/Documents/Experiments/ISMIR_2013/grad/'; % In lunix platform
        dlm    = '/';
end

EXT_TYPE = 2;
switch (EXT_TYPE)
    case 1
        dir = strcat(dir,'pca',dlm);
    case 2
        dir = strcat(dir,'rbm',dlm);
        
        hidNum = [100 500 1000 1200];
        lr_1   = [0.5 0.7];
        lr_2   = [0.7];
        mmt    = [0.1];
        cost   = [0.00002];
    otherwise            
        dir = strcat(dir,'none',dlm);
end

w_num = size(ws,2);

for iiii = 1:200 % set the higher range to search for better features in case of ext using rbm
    log_file = strcat(dir,'exp',num2str(iiii),'.mat')
    inx = resume_from_grid(log_file,8 + w_num);
    if inx(end-w_num+1:end)==ones(1,w_num)
        max_= zeros(1,w_num); 
    else
        max_ = inx(end-w_num+1:end);
    end

    results = zeros(1,w_num);
    W_max = cell(1,w_num);
    vB_max = cell(1,w_num);
    hB_max = cell(1,w_num);
    Ws_max = cell(1,w_num);

for hi = inx(1):size(hidNum,2)
for l1i = inx(2):size(lr_1,2)
% for l1i = inx(3):size(lr_2,2)
for mi  = inx(4):size(mmt,2)
for ci  = inx(5):size(cost,2)    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Feature extraction
features = raw_features;
switch (EXT_TYPE)   
    case 1  % Using PCA        
        assert(~exist('OCTAVE_VERSION'),'This script cannot run in octave');
        coeff = princomp(raw_features);
        coeff = coeff(:,1:6);  % best = 6
        features = raw_features*coeff;
          % normalizing
        mm = minmax(features')';
        inn= (find(mm(1,:)~=mm(2,:)));
        mm = mm(:,inn);
        features = features(:,inn);
        features = (features-repmat(mm(1,:),size(features,1),1))./(repmat(mm(2,:),size(features,1),1)-repmat(mm(1,:),size(features,1),1));
    case 2  % Using rbm      
        conf.hidNum = hidNum(hi);
        conf.eNum   = 100;
        conf.sNum   = size(raw_features,1);
        conf.bNum   = 1;
        conf.gNum   = 1;
        conf.params = [lr_1(l1i) lr_1(l1i) mmt(mi) cost(ci)];
        conf.N    = 50;
        conf.MAX_INC = 10;
        W1 = zeros(0,0);
        [W1 vB1 hB1] = training_rbm_(conf,W1,raw_features);
        features = logistic(raw_features*W1 + repmat(hB1,conf.sNum,1));
    otherwise
        % normalizing
%         mm = minmax(features')';
%         inn= (find(mm(1,:)~=mm(2,:)));
%         mm = mm(:,inn);
%         features = features(:,inn);
%         features = (features-repmat(mm(1,:),size(features,1),1))./(repmat(mm(2,:),size(features,1),1)-repmat(mm(1,:),size(features,1),1));
end

for wi = inx(6):w_num
%% Sub-euclidean computation
w = ws(wi);   % w = subspace window size
num_case = size(trn_inx,1); 
[trnd_12 trnd_13] = subspace_distances(trn_inx,features,indices,w,0);
[tstd_12 tstd_13] = subspace_distances(tst_inx,features,indices,w,0);
cr_ = 0;   % correct rate for training
cr  = 0;   % correct rate for testing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CODE HERE                               %%
[Ws cr_] = gradient_ascent(trnd_12,trnd_13,0.1,0.1,0.00002);

for i = 1:num_case       
 cr = cr + sum((tstd_13{i}-tstd_12{i})*Ws{i}' > 0, 1)/size(tstd_12{i},1);
end
cr = cr/num_case;
if cr_>max_(wi)
    max_(wi) = cr_;
    results(wi)  = cr;
    if EXT_TYPE==2
        W_max{wi} = W1;
        vB_max{wi} = vB1;
        hB_max{wi} = hB1;
        Ws_max{wi} = Ws;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('[window|train|test]= %2d |%f |%f\n',w,cr_,cr);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Using the logging function to save paramters
% and the result for plotting or in grid search
switch EXT_TYPE
    case 1
       % logging(log_file,[100 100 100 100 100 wi cr_ cr max_]);  
    case 2
        logging(log_file,[hi l1i l1i mi ci wi cr_ cr max_ conf.hidNum conf.eNum conf.params]);  
    otherwise
        logging(log_file,[100 100 100 100 100 wi cr_ cr max_]);  
end
end
inx(6)=1;
end
inx(5) = 1;
end
inx(4) = 1;
end
inx(2) = 1;
end
inx(1) = 1;
%% Test on best features

save(strcat(dir,'res_',num2str(iiii),'.mat'),'max_','results','W_max','vB_max','hB_max','Ws_max','ws');
[dummy pos] = max(max_);
fprintf('Accuracy (RBM best fts): w = %d train = %f test = %f\n',ws(pos),max_(pos),results(pos));
clc;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
