%% Function to recover single column of mixing matrix
function [u,peak_val,diff,C,independence_score_m,independence_score_s,zmin,Lmin,Rmin] = ICA_power_with_check(X,maxruns,epsilon, ...
                                                                            B,verbose, ...
                                                                            contrast_fn, ...
                                                                            independence_check, ...
                                                                            independence_check_trials, ...
                                                                            u_init, ...
                                                                            C_init, ...
                                                                            random_projections, ...
                                                                            num_trials, ...
                                                                            M, ...
                                                                            A, ...
                                                                            B1, ...
                                                                            column, ...
                                                                            delta_flag)
% Inputs : X - data with shape (n,k)
%          maxruns - number of iterations to run ICA
%          B - true mixing matrix, used only for debugging purposes
%          verbose - flag to print debug messages
%          contrast_fn - Constrast function to use for ICA, for example
%                        cgf or kurtosis

if nargin==3
    verbose=0;
end

addpath(genpath('../algorithms/'));
addpath(genpath('../util/'));

[~,k] = size(X);

independence_scores = zeros(independence_check_trials, 1);
projections         = random_projections;

if independence_check == 1
    min_independence_score = inf;
    min_independence_score_sd=0;
    itr = 1;
    for r = 1:independence_check_trials
        C_init = contrast_fn.estimate_C(1);
        [u_temp,peak_val_temp,diff_temp,C_temp] = ICA_power(X, maxruns, ...
                                              epsilon, B, verbose, ...
                                              contrast_fn, ...
                                              u_init, ...
                                              C_init, ...
                                              M);
        if(delta_flag)
            Atilde = A; Btilde = B1;
            Atilde(:,column)=u_temp';
            u1=pinv(C_init)*Atilde(:,column);
            Btilde(column,:)=u1/(u1'*Atilde(:,column));
            [independence_score_m,independence_score_s,z,L,R] = Delta_incremental_all(X,Atilde(:,1:column),Btilde(1:column,:),projections,num_trials);
        else
            [independence_score_m,independence_score_s,z,L,R] = Delta_incremental(X,C_temp,u_temp,projections,num_trials);
        end
        
        if(independence_check_trials > 1)
            u_init = randn(k,1);
            u_init = u_init/norm(u_init);
        end
        %independence_score = mean(all_scores);
        %independence_scores(r) = independence_score_m;
        if( independence_score_m+independence_score_s  < min_independence_score-min_independence_score_sd )
            min_independence_score = independence_score_m;
            min_independence_score_sd = independence_score_s;
            zmin = z;
            Lmin = L;
            Rmin = R;
            u = u_temp;
            peak_val = peak_val_temp;
            diff = diff_temp;
            C = C_temp;
        end
        
        if(verbose)
            fprintf("Independence Check Trial - %d, " + ...
                    "Peak Value found - %f, " + ...
                    "Independence Score - %f\n", ...
                    itr, peak_val_temp(end), independence_score);
            keyboard;
        end
        itr = itr + 1;
    end
else 
    [u,peak_val,diff,C] = ICA_power(X, maxruns, epsilon, B, ...
                                    verbose, contrast_fn, u_init, ...
                                    C_init, ...
                                    M);
    if(verbose)
            fprintf("Peak Value found - %f, " + peak_val(end));
    end
end