function [new_w, new_h] =  newhosts_NMF_dist(w,h, newd_out, newd_in, exist_lam, D, landmarks, hosts, weight_on)
% [new_w, new_h] =  newhosts_NMF(w,h, newd_out, newd_in, exist_lam)
% w: Nxd, h: dxN. The position vectors of all landmarks

% N: number of landmarks
% M: number of new hosts
% newd_out: Mxn_landmark out distance
% newd_in: n_landmarkxM in distance

% exist_lam: indexes of measured landmarks
% return: new_w: Mxd, new_h: dxM, such that new_w * h = newd

    error_upbound = 0.5;

    [N,d] = size(w);
    [M,n_landmark] = size(newd_out);
    % nnls, mldivide

    w = w(exist_lam,:);
    h = h(:,exist_lam);    
    
    new_w = zeros(M, d);
    new_h = zeros(d, M);
    
    weight_w = zeros(1, M);
    weight_h = zeros(1, M);
    
    reference_nodes = zeros(M, N);
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
   
    for round = 1:1
        
    for i=1:M        
          
        
        if (round==1)
        rand_existing = randperm(i-1+N);
        dist_vec = zeros(1, i-1+N);
        for j = 1:i-1+N
            if (rand_existing(j) <= N)
                dist_vec(j) = D(hosts(i), landmarks(rand_existing(j)));
            else
                dist_vec(j) = D(hosts(i), hosts(rand_existing(j)-N)) ;
            end
        end
        for j = 1:0
            min_dist = 10000; min_dist_index = 1;
            for k = j+1:(i-1+N)
                if (dist_vec(k) < min_dist)
                    min_dist = dist_vec(k);
                    min_dist_index = k;
                end
            end
            % Swap j and min_dist_index %
            temp = rand_existing(j); rand_existing(j) = rand_existing(min_dist_index); rand_existing(min_dist_index) = temp;
            temp = dist_vec(j); dist_vec(j) = dist_vec(min_dist_index); dist_vec(min_dist_index) = temp;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end
        reference_hosts(i, :) = rand_existing(1:N);
        dist_vec;
        end

        
        % Get the newd_out_tmp, newd_in_tmp, w_tmp, h_tmp
        w_tmp = zeros(N, d);
        h_tmp = zeros(d, N);
        d_out_tmp = zeros(1, N);
        d_in_tmp = zeros(N, 1);
        weight_vec_w = ones(N, 1);
        weight_vec_h = ones(N, 1);
        for j = 1:N
            if (reference_hosts(i,j) <= N)
                w_tmp(j, :) = w(reference_hosts(i,j), :);
                h_tmp(:, j) = h(:, reference_hosts(i,j));
                d_out_tmp(j) = D(hosts(i), landmarks(reference_hosts(i,j)));
                d_in_tmp(j) = D(landmarks(reference_hosts(i,j)), hosts(i));
                if (weight_on == 1)
                    weight_vec_w(j) = 1*rand;
                    weight_vec_h(j) = 1*rand;
%                     weight_vec_w(j) = 1/d_out_tmp(j);
%                     weight_vec_h(j) = 1/d_in_tmp(j);
                end

            else
                w_tmp(j, :) = new_w(reference_hosts(i,j)-N, :);
                h_tmp(:, j) = new_h(:, reference_hosts(i,j)-N);
                d_out_tmp(j) = D(hosts(i), hosts(reference_hosts(i,j)-N));
                d_in_tmp(j) = D(hosts(reference_hosts(i,j)-N), hosts(i));                
                if (weight_on == 1)
                    weight_vec_w(j) = weight_w(reference_hosts(i,j)-N)*rand;
                    weight_vec_h(j) = weight_h(reference_hosts(i,j)-N)*rand;
%                     weight_vec_w(j) = 1/d_out_tmp(j);
%                     weight_vec_h(j) = 1/d_in_tmp(j);
                end
            end
        end
%         min_w = min(weight_vec_w);
%         min_h = min(weight_vec_h);
%         for j=1:N
%             if (weight_vec_w(j) == min_w)
%                 weight_vec_w(j) = eps;
%             end
%             if (weight_vec_h(j) == min_h)
%                 weight_vec_h(j) = eps;
%             end
%         end
        
        
        if (mod(i, 100) == 0)
            fprintf('.');
        end
%         t=h'\newd_out(i,exist_lam)';
%         new_w(i,:)=t';
%         new_h(:,i)=w\newd_in(exist_lam,i);
%         h_tmp
%         d_out_tmp

%          t=h_tmp'\d_out_tmp';
%          new_w(i,:)=t';
%          new_h(:,i)=w_tmp\d_in_tmp;
         
         t = weight_lsq(h_tmp', d_out_tmp', weight_vec_h);
         new_w(i,:)=t';
         new_h(:, i) = weight_lsq(w_tmp, d_in_tmp, weight_vec_w);
         
         % Calculate the error weight %
         weight_w(i)=(error_upbound-median(min(error_upbound, abs(h_tmp'*new_w(i,:)'-d_out_tmp')./(d_out_tmp+0.1)')'))/error_upbound;
         weight_h(i)=(error_upbound-median(min(error_upbound, abs(w_tmp*new_h(:,i)-d_in_tmp)./(d_in_tmp+0.1))'))/error_upbound;
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
    end
    
    fprintf('\n');