function [value,survive,w,it,time] = barycost(F,c,eff,l)
% input arguments:
% F ...  (mxn) matrix of regression vectors, where m is the dimension 
% of the parameter, n is the size of the design space
% c ...  (nx1) vector of costs of design points; 'c' must contain at least
% one value >1 and one value <1, or all the values must be equal to 1
% eff ... required efficiency
% l ...   number of iterations between successive deleting
%      
% output arguments:
% value ... optimal value of det(M)^(1/m)
% survive ... (?x1) vector of indicies of points that have 'survived', i.e. 
% have not been deleted during the process
% w ... (nx1) vector of optimal design weights
% it ... total number of iterations performed
% time ... total time of the computation

tic;
epsilon=1e-06; % used to correct small numerical errors

c(abs(c-1)<epsilon)=1;
szF=size(F); m=szF(1); n=szF(2); nn=n;
% 'n' may vary during the process, 'nn' will keep the constant initial size

% we will re-arrange vector 'c' and matrix 'F' in order to simplify operations
ind=(1:n)'; ind=[ind(c>1);ind(c<1);ind(c==1)];
c=c(ind); F=F(:,ind);

survive=(1:n)';
% numbers of elements in X_+, X_- and X_0, respectively:
n_pl=sum(c>1); n_mi=sum(c<1); n_0=sum(c==1);
n_til=n_pl*n_mi+n_0;  
delta_pl=c(c>1)-1; delta_mi=1-c(c<1);
sums=kron(delta_pl.^(-1),ones(1,n_mi))+kron(ones(n_pl,1),delta_mi.^(-1)');

% initial design
nom_pl=kron(delta_mi',ones(n_pl,1));
nom_mi=kron(delta_pl',ones(n_mi,1)); 
den_pl=kron(delta_pl,ones(1,n_mi)) + kron(delta_mi',ones(n_pl,1));
den_mi=kron(delta_mi,ones(1,n_pl)) + kron(delta_pl',ones(n_mi,1)); 
w_pl=1/n_til*sum(nom_pl./den_pl,2);
w_mi=1/n_til*sum(nom_mi./den_mi,2);
w_0=1/n_til*ones(n_0,1);
w=[w_pl;w_mi;w_0];

M=F*diag(w)*F'; % (mxm) information matrix
d=diag(F'/M *F); % (nx1) variance function
%variance function in points of X_+, X_- and X_0, respectively:
d_pl=d(1:n_pl); d_mi=d((n_pl+1):(n_pl+n_mi)); d_0=d((n_pl+n_mi+1):n); 

a_pl=d_pl./delta_pl; a_mi=d_mi./delta_mi;
% (n_pl x n_mi) matrix of weighted variances:
d_til= (kron(a_pl,ones(1,n_mi))+kron(ones(n_pl,1),a_mi'))./sums; 

MAX=max([d_til(:);d_0]);
it=0; % it-th iteration

while(m/MAX < eff) % main cycle
  it=it+1;
  % computing the new design: 
  wdelta_pl=delta_pl.*w_pl; wdelta_mi=delta_mi.*w_mi;
  S=sum(wdelta_mi);
  w_pl=w_pl.*sum(repmat(wdelta_mi',n_pl,1).*d_til,2)/(m*S);
  w_mi=w_mi.*sum(repmat(wdelta_pl',n_mi,1).*d_til',2)/(m*S);
  w_0=w_0.*d_0/m;
  w=[w_pl;w_mi;w_0];
  
  M=F*diag(w)*F'; 
  d=diag(F'/M *F); 
  d_pl=d(1:n_pl); d_mi=d(n_pl+1:(n_pl+n_mi)); d_0=d((n_pl+n_mi+1):n);
  a_pl=d_pl./delta_pl; a_mi=d_mi./delta_mi;
  d_til= (kron(a_pl,ones(1,n_mi))+kron(ones(n_pl,1),a_mi'))./sums; 
  MAX=max([d_til(:);d_0]);
  
  %deleting:
  if (mod(it,l) == 0)
     eps=MAX-m;
     h_m=m*(1+eps/2-sqrt(eps*(4+eps-4/m))/2);
     
     % logical vectors implying which points of X_+, X_-, X_0 are deleted:
     del_pl=logical(max(d_til,[],2)<h_m-epsilon); %(n_plx1)
     del_mi=logical(max(d_til,[],1)'<h_m-epsilon); %(n_mix1)
     del_0=logical(d_0<h_m-epsilon); %(n_0x1)  
     del=[del_pl;del_mi;del_0];  
                
     if (sum(del)>0) % if at least one point is deleted
        survive(del)=[];
        F(:,del)=[];
        c(del)=[];
        n=n-sum(del); 
        n_pl=n_pl-sum(del_pl); n_mi=n_mi-sum(del_mi);  n_0=n_0-sum(del_0);     
        w_pl(del_pl')=[]; w_mi(del_mi')=[];  w_0(del_0')=[]; 
        delta_pl=c(c>1)-1; delta_mi=1-c(c<1); 
        sums(del_pl,:)=[]; sums(:, del_mi)=[];
          
        % re-normalizing the design:
        w_pl=w_pl(:); w_mi=w_mi(:); w_0=w_0(:);
        
        s_pl=sum(w_pl); s_mi=sum(w_mi); s_0=sum(w_0); s=s_pl+s_mi+s_0;
        sd_pl=sum(delta_pl.*w_pl); sd_mi=sum(delta_mi.*w_mi);
           
        h_pl=sd_mi*(s_pl + s_mi)/(s*(s_pl*sd_mi + s_mi*sd_pl));
        h_mi=sd_pl*(s_pl + s_mi)/(s*(s_pl*sd_mi + s_mi*sd_pl));
        h_0=1/s;
        w_pl=h_pl*w_pl; w_mi=h_mi*w_mi; w_0=h_0*w_0;
        w=[w_pl;w_mi;w_0];
          
        % updating:        
        M=F*diag(w)*F'; 
        d=diag(F'/M *F); 
        d_pl=d(1:n_pl); d_mi=d((n_pl+1):(n_pl+n_mi)); d_0=d((n_pl+n_mi+1):n);
        a_pl=d_pl./delta_pl; a_mi=d_mi./delta_mi;
        d_til= (kron(a_pl,ones(1,n_mi))+kron(ones(n_pl,1),a_mi'))./sums; 
        MAX=max([d_til(:);d_0]);
     end 
  end 
end

value=det(M)^(1/m);
survive=ind(survive);
w2=w; w=zeros(nn,1); w(survive)=w2;
survive=sort(survive);
time=toc;

end