function [out_host, in_host] = phoneix_simple(D, dim, N, K)


%(D_change, dim, N, K, iteration, Select_neighbor, Select_kernel, Select_stat_neg, Select_delta)

D_change = D;

% Parameters
%C = 5;
round_bound = 30;


display('Phoenix(simple)');
neighbor = zeros(N,K);
error_in = zeros(1,N);
error_out = zeros(1,N);
for i=1:N
    out_host(i,:)=rand(1,dim);
    in_host(:,i)=rand(dim,1);

    % Randomly selected Neighbors
      tmp = randperm(N);
      point = 1;
      for j=1:K
          if (i ~= tmp(point))
              neighbor(i, j) = tmp(point);
              point = point + 1;
          else
             while(D_change(i, tmp(point))== 0 || D_change(i, tmp(point)) < 0)
                point = point + 1;
             end
              neighbor(i, j) = tmp(point);
              point = point+1;
          end
      end
%      neighbor(i,:)= tmp(1:K);
      
      % remove the zero distance...
    
%     tmp = randperm(N-closest);
%     tmp = tmp + closest;
%     
%     Latency_vector = D_change(i, :); % Latency to all nodes from node i
%     [Sorted_Latency_vector, latency_index] = sort(Latency_vector);
%     neighbor(i, 1:closest) = latency_index(1:closest);
%     neighbor(i, closest+1:K) = latency_index(tmp(1:K-closest));

    error_out(i) = 1;
    error_in(i) = 1;
end

% Out : X in the paper, In : Y in the paprt
%delta_out = delta;
%delta_in = delta;

%for j=1:15

w=zeros(N, K, dim);
h=zeros(N, dim, K);
D_host2landmark=zeros(N, K);
D_host2landmark_out=zeros(N, K);
D_host2landmark_in=zeros(N, K);

%for j = 1:100
%for j = 1:(iteration/2)
fpre = [];
npre = [];
for round = 1:round_bound

    if (round == 1)
        % First round : Joining %
        
        % NMF for first K nodes %
        D_landmark = D(1:K, 1:K);
        [out_landmark, in_landmark] = NMF(D_landmark, dim);
        %%%%%%%%%%%%%%%%%%%%%%%%%
        out_host(1:K, :) = out_landmark;
        in_host(:, 1:K) = in_landmark;
        % Setup neighbors %
        for i=1:N
            if (i<=K)
                reference_hosts(i, :) = 1:K;
            else
                rand_existing = randperm(i-1);
                reference_hosts(i, :) = rand_existing(1:K);
            end
        end     
    end
    for i = 1:N
        
        % K: number of landmarks
        % w: Kxd, h: dxK. The position vectors of all landmarks

%        tmp = randperm(N); 
%        neighbor(i,:)= tmp(1:K);
        for j=1:K
            if ((neighbor(i, j) < 0 | D_change(i, neighbor(i, j)) < 0) & round > 1)
                while(1)
                    ok = 1;
                    nb_candidate = ceil(rand*N);
                    if (D_change(i, nb_candidate) < 0)
                        continue;
                    end
                    for jj=1:K
                        if (nb_candidate == neighbor(i, jj))
                              ok = 0;
                              break;
                        end
                    end
                    if (ok == 1)
                          neighbor(i,j)=nb_candidate;
                        break;
                    end
                end
            end
        end
        if (round == 1 & i < K)
            continue; % Already has NC
        end
        temp_w = out_host(neighbor(i,:),:);
        temp_h = in_host(:, neighbor(i,:));
        temp_D_host2landmark = D_change(i, neighbor(i, :));
         temp_D_host2landmark_out = temp_D_host2landmark(1:K);
         temp_D_host2landmark_in = temp_D_host2landmark(1:K);

        if (round == 1) % Initial NC Calculation
            for ii=1:K
                if (D(i,  neighbor(i,ii)) < 0)
                    weight_out_vec(ii) = eps;             
                    weight_in_vec(ii) = eps;           
                else
                    weight_out_vec(ii) = 1;             
                    weight_in_vec(ii) = 1;           
                end
            end
            t = weight_lsqnonneg(temp_h', temp_D_host2landmark_in', sqrt(weight_in_vec)');

            out_host(i, :) = t';
        
            in_host(:, i) = weight_lsqnonneg(temp_w, temp_D_host2landmark_out', sqrt(weight_out_vec)');
        end
             

        


        for ii=1:K
                weight_out_vec(ii) = 1;
                weight_in_vec(ii) = 1;
        end

        
%         temp_w = temp_w(index_nblist_out, :);
%         temp_h = temp_h(:, index_nblist_in);
%         temp_D_host2landmark = temp_D_host2landmark(index_nblist);
         
         temp_D_host2landmark_out = temp_D_host2landmark(1:K);
         temp_D_host2landmark_in = temp_D_host2landmark(1:K);
         
      
        
        % x = lsqnonneg(C,d) 
        % returns the vector x that minimizes norm(C*x-d) subject to x >= 0. C and d must be real.

        %  out_host(i, :) [1*d] * In_NCs[d*K] = D_host2landmark[1*K];        
        %  => out_host(i, :) * h = D_hosts2landmark[1*K];
        %  => h' * out_host(i, :)' = D_host2landmark';
%        size(temp_h)
%        size(temp_D_host2landmark)

%        t = lsqnonneg(temp_h', temp_D_host2landmark_in');
        t = weight_lsqnonneg(temp_h', temp_D_host2landmark_in', sqrt(weight_in_vec)');

        out_host(i, :) = t';
        
        %  Out_NCs[K*d] * in_host(:, i) [d*1] = D_host2landmark'[K*1];
        %  => w * in_host(:, i) = D_hosts2landmark'[K*1];
        
%        in_host(:, i) = lsqnonneg(temp_w, temp_D_host2landmark_out');
        in_host(:, i) = weight_lsqnonneg(temp_w, temp_D_host2landmark_out', sqrt(weight_out_vec)');

%         w = backup_w;
%         h = backup_h;
%         D_host2landmark = backup_D_host2landmark;

         % Randomly kick out one %
%          start_point = ceil(rand * K);
%          tmp_index = start_point;
%          for ii = 1:K
%              if (score_out_vec(tmp_index) == out_max || score_in_vec(tmp_index) == in_max)
%                  neighbor(i, tmp_index) = -1;
%                  break;
%              end
%              tmp_index = tmp_index + 1;
%              if (tmp_index > K)
%                  tmp_index = 1;
%              end
%          end
     
    end;
    
end;

