function [Curve, Strip, Rank, pvals] = ...
    run_multi_line_2D_bioinfo06(sample, asymmetry, num_shfts, min_num, max_level, thresh_global, sc_rat, ratio_beta, toon, solve, do_svd)

% This function applies the MSC algorithm both for cDNA arrays (asymmetry = 0) 
% and to ChIP-on-chip data (asymmetry = 1).
% It follows ideas of the paper "Multiscale Curve and Strip Constructions",
% Lerman, McQuown, Mishra, 2006.
% Original code has been written by Gilad Lerman and Joseph McQuown and it
% has been revised carefully by Guangliang Chen.

% inputs    sample          - real Nx2 matrix of data. It is assumed that
%                             the data is given in MA coordiantes. That is,
%                             one has already applied the transformation:
%                             sample = [log2(sample(:,1)) log2(sample(:,2))];
%                             sample = [0.5*(sample(:,1)+sample(:,2)) -sample(:,1)+sample(:,2)];
%           asymmetry       - integer 0/1. =0 for symmetric data; =1 for
%                             asymmetric data
%           num_shifts      - integer. number of shifts to smooth the curve
%                             and strip ($n_sh$ in LMM06; we recommend num_shfts = 30)
%           min_num         - integer; min number of points with which to fit local
%                             stats (Constant $n_0$ in LMM06). We recommend
%                             min_num = 30.
%           max_level       - integer; max number of levels algorithm is to run
%                             (constant $l_0$ in LMM06). We recommend
%                             max_level = 10.
%           thresh_global   - real; threshold for main stopping time condition.
%                             It is an upper bound for precentage of
%                             outliers (constant $\alpha_0$ in LMM06). See
%                             LMM06 for discussion of its choice.
%           sc_rat          - real; % in cDNA arrays (asymmetry = 0) we use sc_rat to form $c_0$ of LMM06.
%                             $c_0$ is sc_rat times the ratio of the smallest height of a rectangle 
%                             containing the data over the length of the base of that rectangle. If the
%                             sc_rat is 1, then the expression above (=the ratio) is our advocated choice for the parameter $c_0$
%                             In ChIP-on-chip data (asymmetry = 1), we use
%                             it to form $c_{-}$ of LMM06. $c_{-}$ is the
%                             multiplication of $sc_rat$ by the ratio of twice the minimal y-vlaue (height of a triangle 
%                             in a very symmetric case) over the
%                             sidelength of smallest interval containing the $x$-values. If sc_rat=1, then the 
%                             expression above (=the ratio) is our advocated choice for the parameter $c_{-}$.
%           ratio_beta      - real and positive or -1; It dectates the parameter $delta_0$ of LMM06.
%                             If asymmetry = 0, it is the ratio of $\delta_0$/$c_0$.
%                             If asymmetry = 1, it is the ratio of
%                             $\delta_0$/$c_{-}$.
%                             If beta_thresh = -1, then the stopping time
%                             condition involving the parameter $delta_0$
%                             is deactivating. We recommend on that option.
%           toon            - integer 0/1. if =1, then a cartoon shows the process
%           solve           - integer 0/1. =1 if the output strip is
%                             \hat{S} (see LMM06), 0 if it is \tilde{S} (LMM06)
%                             The output is more precise if solve = 1, but
%                             slower.
%          do_svd           - integer 0/1. If do_svd = 0, then no change to
%                             the data (we compute at this part the main regression
%                             line, instead of the following parts, and subtract its value from 
%                             the y-values). If do_svd=1, then shift and rotate it
%                             onto main principal axis.
%
% outputs   Strip           - N*1 real vector. Estimates local conditional standard deviation.
%                             \tilde{S} if solve = 0 and \hat{S} if solve = 1.
%           Curve           - N*(D-1) real vector. The regression curve for
%                             the stable distribution (denoted by C in LMM06).
%           Rank            - N*1 real vector of ranks (see R in LMM06). It
%                             uses \tilde{S} if solve = 0 and \hat{S}
%                             otherwise.
%           pvals           - N*1 real vector of p-values dictated by rank.
%

N = size(sample,1);
order_apprxmtn = 1; %degree of approximation (=1 for local lines)

if do_svd
    [sample_min,idx_min] = min(sample(:,2));
    mean_sample = mean(sample);
    sample = sample - ones(N,1) * mean_sample;
    [u,s,v] = svd(sample,0);
    sample = sample * v;

    if sample(idx_min(1),2) > 0
        sample(:,2) = -sample(:,2);
    end
    %sample = sample + ones(N,1) * mean_sample;
    clear u s sample_min
else
    %[b, err] = get_line_D_svd_bioinfo06(sample,order_apprxmtn);
    b = regress(sample(:,2),[ones(N,1) sample(:,1)]);    
    sample(:,2) = sample(:,2) - [ones(N,1) sample(:,1)]*b;
    %clear b
end

% sorting x-coordinates
[sample_sort(:,1), idx_sort] = sort(sample(:,1)); 
sample_sort(:,2) = sample(idx_sort,2);
inv_idx_sort(idx_sort) = 1:N;

%%%%%%%%%%%%%  Set parameters for the multiscale algorithm %%%%%%%%%%%%%%%%

min_s = sample_sort(1,1);
max_s = sample_sort(N,1);
interval_length = max_s - min_s;

% The following three lines ensure that we do not loose points 
% of sample at the endpoints of the bounding interval.
minC = min_s - 1e-5 * interval_length / 2;  %left boundary of shortest bounding interval
interval_length = (1+1e-5) * interval_length;
maxC = minC + interval_length; %right boundary of shortest bounding interval

% num_shfts = 30;


if num_shfts == 1
   shft = -1;
else
   shft = (0:num_shfts-1) / num_shfts;
   interval_length = 2*interval_length;
end

% if asymmetry = 0, then sc_up and sc_down are both c_0 in LMM06.
% if asymmetry = 1, then sc_up is c_{+} and sc_down is c_{-}
% Remark: We have a factor 2 in defining sc_up and sc_down, because they 
% are The comparison has to be done indeed in the first level,
% the zeroth level of shifts is a bit arbitrary.
% asymmetry = 1;
if asymmetry == 0 %symmetric data
    sc_down = 2*sc_rat*max(abs(sample(:,2)))/interval_length;
    sc_up = sc_down;
else %asymmetric data like CHIP-on-chip data
    sc_up = 2^(max_level+1)*(max(sample(:,2))/interval_length);
    sc_down = 2*sc_rat*abs(min(sample(:,2)))/interval_length;
end


%%%coefficient >1 or <1???

if isempty(ratio_beta) || ratio_beta == -1
    ratio_beta = max(sc_up,sc_down)/sc_down;
end

% We set the value of this parameter to be sufficiently large, so that it
% is not taking into account (the corresponding stopping time aren't
% considered).
beta_thresh = sc_down*ratio_beta;



% The parameter below corrects edge effects (at the boundary interval). 
% One may need to decrease it to 0.5 if getting warning: "dividing by zero"
edge_ratio = 0.8;% controls edge effects


curve_option = 1; % =1 for unnormalized data

partition = 0; %0 partition using midpoint (our method)
               %1 using median of Q_til, 2 using median of Q_hat, 3 for another normalization
               
               
% fprintf(1,'\nmin_num = %6.2f \n', min_num)
% fprintf(1,'max_level = %6.2f \n', max_level)
% fprintf(1,'thresh_global = %6.5f \n', thresh_global)
% fprintf(1,'num_shfts = %6.2f \n', num_shfts)
% fprintf(1,'sc_up = %6.5f \n', sc_up)
% fprintf(1,'sc_down = %6.5f \n', sc_down)
% fprintf(1,'beta_thresh = %6.4f \n', beta_thresh)
% fprintf(1,'edge_ratio = %6.2f \n', edge_ratio)
% fprintf(1,'order_apprxmtn = %6.2f \n', order_apprxmtn)
% fprintf(1,'curve_option = %6.2f \n', curve_option)

% %%% initializing bounding interval

Strip_avg = zeros(N,1); 
Curve_avg = zeros(N,1);  
sum_Edge_weight = zeros(N,1); 
    
%%%%%%%%%%%%%%   Apply repetedly the algorithm with different shifts
for i = 1 : num_shfts
    
    [Strip,Curve,Edge_weight] = ...
        find_multi_line_2D_bioinfo06(sample_sort,minC,maxC,max_level,thresh_global,sc_up,sc_down,...
            min_num,beta_thresh,shft(i),edge_ratio,order_apprxmtn,curve_option,toon,partition,asymmetry,solve);

    Strip_avg = Strip_avg + Edge_weight.*Strip;
    Curve_avg = Curve_avg + Edge_weight.*Curve;
    sum_Edge_weight = sum_Edge_weight + Edge_weight; 

end

clear sample_sort

%may get warning "dividing by zero" due to the following command, can try
%decreasing edge_ratio or increasing sc_down
%if min(sum_Edge_weight) > 0
    Strip_avg = Strip_avg ./ sum_Edge_weight;
    Curve_avg = Curve_avg ./ sum_Edge_weight;
% else
%     Curve_avg = zeros(N,1);
%     Strip_avg = ones(N,1)*sqrt(sum(sample_sort(:,2).^2)/(N-1));
% end

Curve = Curve_avg(inv_idx_sort,:);
Strip = Strip_avg(inv_idx_sort,:);

clear Curve_avg Strip_avg

if asymmetry == 1
    Rank = -(sample(:,2) - Curve)./Strip;
    Rank = max(Rank,0);
else
    Rank = abs(sample(:,2) - Curve)./Strip;
end
%Rank = min(Rank,0);

%One sided pvalues determined here normcdf = 0.5 * erfc(-(x-mu)./(sqrt(2)*ma));
%where erfc(x) = 1 - erf(x)
%p = normcdf(sample(:,2),Curve,2);

%not for chip-on-chip
if asymmetry == 0
    %pvals = normcdf(Rank);
    pvals  = erfc(Rank/sqrt(2));
else
    %pvals = 1-normcdf(Rank);
    pvals  = erfc(Rank/sqrt(2))/2; %CHIP-2-chip
    %pvals  = erfc(Rank/sqrt(2));
end

if ~do_svd
    Curve = Curve + [ones(N,1) sample(:,1)]*b;
end

