function root_index = chooseRoot(data_discrete,noisy_magnitude,m,threshold)

% % % Copyright (c) 2015, Mark Styczynski and Weiwei Yin
% % % All rights reserved.
% % % 
% % % Redistribution and use in source and binary forms, with or without
% % % modification, are permitted provided that the following conditions are
% % % met:
% % % 
% % %     * Redistributions of source code must retain the above copyright
% % %       notice, this list of conditions and the following disclaimer.
% % %     * Redistributions in binary form must reproduce the above copyright
% % %       notice, this list of conditions and the following disclaimer in
% % %       the documentation and/or other materials provided with the distribution
% % % 
% % % THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% % % AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% % % IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% % % ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% % % LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% % % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% % % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% % % INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% % % CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% % % ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% % % POSSIBILITY OF SUCH DAMAGE.


[numOfobs, numOfvar] = size(data_discrete) ; 

%create noisy data 
for ll = 1:numOfvar
    for nn = 1:noisy_magnitude

        noisy_training(:, noisy_magnitude*(ll-1)+nn) = data_discrete(randperm(numOfobs), ll) ; 
    end 
end 


%combine real and noisy data
training_whole = [data_discrete noisy_training] ; 

% calculate s-distribution
for ll = 1:numOfvar

    ll;

    for kk = 1:(numOfvar*(noisy_magnitude+1))

        kk;

        [s(kk), flag(kk) ] = SignificanceTestByPermutation(training_whole(:, ll), training_whole(:, kk), m, threshold); 

    end 

    % calculate the p-values by comparing the s-distributions for real
    % and noisy data
    s_record(ll, :) = s ; 

    s(ll) = [] ; % remove the s for the root itself. 

    [h(ll), p(ll)] = ttest2(s(1:(numOfvar-1)), s(numOfvar:end), 0.05, 'both', 'unequal') ; % test with original s-values

%         [h(ii, jj), p(ii, jj)] = ttest2(abs(temp(1:36)), abs(temp(37:end)), 0.05, 'both', 'unequal') ; % test with the absolute values of  s

end 

% pick the node with the smallest p-value as the root
[min_value_p, root_index] = min(p) ;

% after record the selected root in root_index, 
s_record(root_index, root_index) = 0  ; 



