function [cost_SID_CSE, CSE_table, key, p, shift, porm, CSE_bp] = SID_CSE(ch,wl, maxncse)

%SID_DMST + CSE (Common Subexpression Elimination)
%Yongtao Wang, yongtaow@gmail.com
%
%For info about SID_DMST, refer to the help info of function SID_DMST.m
%
%INPUTS:
%   ch: filter coefficients. For linear-phase filter, contains only the unique coefficients
%   wl: wordlength
%   maxncse: the maximum number of common subpexressions that can be shared
%     if maxncse=-1 means that there is no upper bound
%   logicdepth: upper limit on the logic depth of the final implementation
%OUTPUTS:
%   cost_SID_CSE: number of adders required for the multiplier block of the filter
%   CSE_table: a sructural array, each element containing a subexpression in field 'bp'
%     CSE_table(1).bp is always 1, CSE_table(2).bp is represented by '2', CSE_table(3).bp by 3 and so on. 
%   p: parent of each node from 1 to M. p(M)=0 since it is the root of DMST
%   key: weight of the edges in the DMST. key(i) is weight of the edge directed from p(i) to node i
%   shift: the amount of shift if a SID coefficient is used, see MRPF_DMST for more explanation
%   porm: "Plus OR Minus" in forming the SID coeff., see MRPF_DMST for more explanation
%   CSE_bp: from p,shift and porm we can determine the coefficients used. For example, for j-th coefficient in ch,
%      dc_j = ch(j)+porm(j)*2^shift(j)*ch(p(j)) is used for implementing ch(j)*x, where x is the input of the filter.
%      So CSE_bp(j,:) is the bit pattern of dc_j after CSE.
%      For example, dc_j=[0 1 0 1 0 0 0 -1 0 1] before CSE. Assume CSE_table(2).bp=[1 0 1] and CSE_table(3).bp=[1 0 -1],
%      then dc_j=[0 2 0 0 0 0 0 -3 0 0].
%KEY CAREABOUTS:
%1. From p, shift, porm, CSE_table, CSE_bp, we have all necessary info for
%filter implementation.
%
%EXAMPLE:
% F = [0.4 .8]; A = [1 0]; wp=0.005; ws=0.001;
% DEV = [wp ws]; [n, fo, mo, w] = remezord(F,A,DEV);
% c = remez(n, fo, mo, w);
% wl = 16;
% cq=[];
% for i=1:length(c)
%     cq(i) = quantize(c(i), wl, 1);
% end
% ch = cq(1:ceil(length(c)/2)); 
% %number of adders with simple CSD implementation
% cost_CSD=0; 
% for i=1:length(ch)
%     cs = CSD(ch(i),wl);
%     a = sum(abs(cs));
%     if a>0
%         cost_CSD = cost_CSD + a -1;
%     end
% end
% maxncse=-1;
% [cost_SID_CSE, CSE_table, key, p, shift, porm, CSE_bp]= SID_CSE(ch,wl, maxncse);

%construction of the directed multiple graph
M = length(ch)+1;
G = cell(M, M); GM={M,M};
shift=[]; porm=[];
%G{i,j}: edge directed from cj to ci and SID coeefficients ci+/-2^L*cj
%edges from cM to ci (i in [1, M-1]): original coefficients
for i=M
    for j=1:M-1
        edge = ch(j);  
        cbstr = CSD(edge, wl); %change CSD so that it outputs a vector consisting of 0, 1, and -1
        cost = sum(abs(cbstr));
        shift=0; %record the value of shift 
        porm=0; %plus or minus: 0 means that there is no differential coeff.
        if cost > 0
            cost = cost - 1;
        end
        G{j,i} = [edge; cost; shift; porm];
        ebpm = cbstr; %edge bit-pattern matrix
        GM{j,i} = ebpm;
    end
end
for i=1:M-1
    chb = CSD(ch(i), wl); ht = htdetect(chb);
    for j=1:M-1
        if i~=j
            edge = []; cost = []; shift=[]; porm=[];
            for L = 0 : ht(1)-1
                if abs(ch(j)-2^L*ch(i)) < 4/3
                    edge = [edge ch(j)-2^L*ch(i)];
                    porm=[porm -1];
                    shift = [shift L];
                end
                if abs(ch(j)+2^L*ch(i)) < 4/3
                    edge = [edge ch(j)+2^L*ch(i)];
                    porm=[porm 1];
                    shift = [shift L];
                end
            end
            for L = 1 : wl-ht(2)
                if abs(ch(j)-2^(-L)*ch(i)) < 4/3
                    edge = [edge ch(j)-2^(-L)*ch(i)];
                    porm=[porm -1];
                    shift = [shift -L];
                end
                if abs(ch(j)+2^(-L)*ch(i)) < 4/3
                    edge = [edge ch(j)+2^(-L)*ch(i)];
                    porm=[porm 1];
                    shift = [shift -L];
                end 
               
            end
            ebpm=[]; %edge bit pattern matrix
            for k = 1 : length(edge)
                cbstr = CSD(edge(k), wl);
                cost(k) = sum(abs(cbstr));
                ebpm(k,:) = cbstr;
            end
            GM{j,i} = ebpm;
            G{j,i} = [edge; cost; shift; porm]; 
        end
    end
end
%non-existent edges defined as edge=-100 and cost=1000
for i=1:M
    edge=[]; cost=[];
    edge = -100; cost = 1000; 
    shift=-100; porm=0;
    GM{i,i}=[];
    G{i,i} = [edge; cost;shift;porm];
end
for i=1:M-1
    edge=[]; cost=[];
    edge = -100; cost = 1000;
    shift=-100; porm=0;
    GM{M,i}=[];
    G{M,i} = [edge; cost;shift;porm];
end
%note: G_cost is a MxM matrix
%selected_edge_index: index of the selected edge 
G_cost=[]; selected_edge_index=[];
for i=1:M
    for j=1:M
        [Y, I] = min(G{i,j}(2,:));
        G_cost(i,j)=Y;
        selected_edge_index(i,j) = I(1);
    end
end
%DMST algorithm: cost without doing CSE
[key, p] = DMST(G_cost);

%Apply CSE: steep descent algorithm + recursive elimination of 2-bit subexpressions

BPT = cell(M,M); %store the detailed BP info. for all the edges of the graph
%each element of BPT is a structured array MEBPT storing the BP info about the multiple edges joinging i and j
CSE_table=[];
noCSE = 0; cti=1; CSE_table(1).bp = [1];
while ~noCSE
    pp=1; GBPT_tmp = []; 
    for i=1:M
        for j=1:M
            MEGBPT=[]; q=1; MEBPT=[];
            for k=1:size(GM{i,j},1) %multiple edges joining vertices i and j
                cb = GM{i,j}(k,:);
                MEBPT(k).LBPT = bpdetect_Har(cb);
                %only consider edges whose current cost is no greater than G_cost(i,j)+1 
                %since other edges can not be part of the solution
                if nzbitcount(cb) <= G_cost(i,j)+1 
                    for m = 1:length(MEBPT(k).LBPT)
                        [IsIn, index] = IsInBPT_Har(MEBPT(k).LBPT(m).bp, MEGBPT);
                        if IsIn
                            %for a bp occurs in multiple edges, use the maximum frequency, not count all occurences
                            MEGBPT(index).freq = max(MEGBPT(index).freq, MEBPT(k).LBPT(m).freq); 
                        else
                            MEGBPT(q).bp = MEBPT(k).LBPT(m).bp; MEGBPT(q).freq = MEBPT(k).LBPT(m).freq;
                            q=q+1;
                        end
                    end 
                    for m = 1: length(MEGBPT)
                        [IsIn, index] = IsInBPT_Har(MEGBPT(m).bp, GBPT_tmp);
                        if IsIn
                            GBPT_tmp(index).freq = GBPT_tmp(index).freq + MEGBPT(m).freq; 
                        else
                            GBPT_tmp(pp).bp = MEGBPT(m).bp; GBPT_tmp(pp).freq = MEGBPT(m).freq;
                            pp=pp+1;
                        end
                    end
                end
            end
            BPT{i,j} = MEBPT;
        end
    end
    pp=1; freq=[]; GBPT=[];
    for i =1:length(GBPT_tmp)
        if GBPT_tmp(i).freq >=2
            GBPT(pp).bp = GBPT_tmp(i).bp;          
            freq=[freq GBPT_tmp(i).freq]; 
            pp=pp+1;
        end
    end
    if length(freq)==0
        noCSE = 1;
    else 
        cost_CSE=[];
        for bpindex =1:length(freq)
            bp = GBPT(bpindex).bp;
            cti = cti + 1; 
            %update the edge weights assuming bp is eliminated
            GM_tmp=GM; G_tmp=G;
            for i=1:M
                for j=1:M
                    MEBPT = BPT{i,j}; %multiple edge bit pattern table
                    for k=1:length(MEBPT)
                        cb = GM{i,j}(k,:);
                        [IsIn, index] = IsInBPT_Har(bp, MEBPT(k).LBPT);
                        if IsIn
                            G_tmp{i,j}(2,k) = G_tmp{i,j}(2,k)- size(MEBPT(k).LBPT(index).bpi,1);
                        end
                    end
                end
            end
            G_cost_tmp=zeros(M); %MxM matrix
            for i=1:M
                for j=1:M
                    G_cost_tmp(i,j)=min(G_tmp{i,j}(2,:));
                end
            end
            %DMST algorithm
            [key_tmp, p_tmp] = DMST(G_cost_tmp);
            cost_CSE(bpindex) = sum(key_tmp)+1; cti = cti-1;
        end
        %steep descent
        [Y, I] = min(cost_CSE);
        if Y>=sum(key) %no reduction: stop the process by setting noCSE=1
            noCSE=1;
        else
            bp = GBPT(I(1)).bp; %steep descent
            cti = cti + 1;
            %really eliminate the selected bp
            GM_tmp=GM; G_tmp=G;
            for i=1:M
                for j=1:M
                    MEBPT = BPT{i,j};
                    for k=1:length(MEBPT)
                        cb = GM{i,j}(k,:);
                        if nzbitcount(cb) <= G_cost(i,j)+1 
                            [IsIn, index] = IsInBPT_Har(bp, MEBPT(k).LBPT);
                            if IsIn
                                for t = 1: size(MEBPT(k).LBPT(index).bpi,1)
                                    if bp(1) == GM_tmp{i,j}(k,MEBPT(k).LBPT(index).bpi(t, 1)) % exactly same
                                        GM_tmp{i,j}(k,MEBPT(k).LBPT(index).bpi(t, 1))= cti;
                                    else %same as the complementary bp
                                        GM_tmp{i,j}(k,MEBPT(k).LBPT(index).bpi(t, 1))= -cti;
                                    end
                                    for h = 2:size(MEBPT(k).LBPT(index).bpi, 2)
                                        GM_tmp{i,j}(k,MEBPT(k).LBPT(index).bpi(t, h))=0;
                                    end
                                    G_tmp{i,j}(2,k) = G_tmp{i,j}(2,k)-1;
                                end
                            end
                        end
                    end
                end
            end
            G_cost=[];selected_edge_index_tmp=[];
            for i=1:M
                for j=1:M
                    [YG, IG] = min(G_tmp{i,j}(2,:));
                    G_cost(i,j)=YG;
                    selected_edge_index_tmp(i,j) = IG(1);
                end
            end
            %DMST algorithm
            [key_tmp, p_tmp] = DMST(G_cost); 
            CSE_table(cti).bp = bp;
            GM = GM_tmp; G = G_tmp;
            key = key_tmp; p = p_tmp; 
            selected_edge_index = selected_edge_index_tmp;
        end
    end
end

shift=[]; porm=[];
CSE_bp=[]; %bit pattern of selected edge with subexpression eliminated
for i=1:length(p)-1
    shift(i) = G{i,p(i)}(3,selected_edge_index(i,p(i)));
    porm(i) = G{i,p(i)}(4,selected_edge_index(i,p(i)));
    CSE_bp(i,:) = GM{i,p(i)}(selected_edge_index(i,p(i)),:);
end
cost_SID_CSE = sum(key)+length(CSE_table)-1;
