function [fwdiff]=COMMU_fdiff(dataX,dataY,same_index,diff_index,Dl,w,B,C)
 %%compute the minimizer in Eq.9 in the article
[n,p]=size(dataX);
Ks=size(same_index,2);
Kr=size(diff_index,2);
q=size(dataY,2);
Dm=eye(p,p);
for i=1:q
 Dm=Dm+w(i)*B(:,i)*B(:,i)';%p*p
end 

temp=zeros(1,q);
D=zeros(n,Kr*Ks); 
for k=1:Kr    
     Dr=dataX-dataX(diff_index(:,k),:);
     for j=1:Ks
          Ds=dataX-dataX(same_index(:,j),:);   
          %Dld=dataY(same_index(:,j),:)-dataY(diff_index(:,k),:);
          D(:,Ks*(k-1)+j)=diag(dataY(same_index(:,j),:)*Dl*dataY(diff_index(:,k),:)')-diag(Dr*Dm*Dr')+diag(Ds*Dm*Ds');               
     end
end
[maxD,index]=max(D,[],2);%Ks*(k-1)+j=index

k_r=floor(index/Ks)+1;

%kr=(index/Ks==fix(index/Ks))?(index/Ks):(floor(index/Ks)+1);
k_s=mod(index,Ks);
k_s(k_s==0)=Ks;
k_r(k_s==Ks)=k_r(k_s==Ks)-1;
big_zero=find(maxD>0);

for i=1:length(big_zero)    
   for m=1:q                
     Btemp=B(:,m)*B(:,m)';
     temp(:,m)=temp(:,m)+maxD(big_zero(i))*((dataX(big_zero(i),:)-dataX(same_index(big_zero(i),k_s(big_zero(i))),:))*Btemp*(dataX(big_zero(i),:)-dataX(same_index(big_zero(i),k_s(big_zero(i))),:))'-(dataX(big_zero(i),:)-dataX(diff_index(big_zero(i),k_r(big_zero(i))),:))*Btemp*(dataX(big_zero(i),:)-dataX(diff_index(big_zero(i),k_r(big_zero(i))),:))');
    %temp(:,m)=temp(:,m)+diag(maxD)*((dataX-dataX(same_index(:,k_s),:))*Btemp*(dataX-dataX(same_index(:,k_s),:))'-(dataX-dataX(diff_index(:,k_r),:))*Btemp*(dataX-dataX(diff_index(:,k_r),:))');
   end
end
  
   fwdiff=2*C/n*temp;
end