function [fitted, coordone,coordtwo,addcon,fval,vaf,exitflag] = ...
    lpfit(prox,permone,permtwo,start_vector,p)

%lpfit.m provides the two-dimensional coordinates used to construct lp
%distances fit (in a confirmatory manner) to the given proximities.

%The inputs:
%
%the n x n symmetric dissimilarity matrix with zero main diagonal entries;
%two fixed permutations,  PERMONE and PERMTWO, of the objects along the
%two dimensions;
% the value of p to obtain the l_{p} distances (p = 1: city-block distances; p
% = 2: Euclidean distances; large p: dominance distances).
%the START_VECTOR of length 2n+1.

%The outputs:
%
%the values fitted to the proximity matrix are in FITTED;
%the object coordinates along the two dimensions, COORDONE and COORDTWO;
%the additive constant, ADDCON;
%the values of the least-squares criterion, FVAL, in finding the
%coordinates;
%the variance-accounted-for measure, VAF;
%the exit condition, EXITFLAG, from using the MATLAB optimization routine,
%fmincon.m.
%
% syntax: [fitted, coordone,coordtwo,addcon,fval,vaf,exitflag] = ...
%    lpfit(prox,permone,permtwo,start_vector,p)




n = size(prox,1);

coordone = zeros(n,1);
coordtwo = zeros(n,1);

A = zeros(2*(n-1),2*n+1);
b = zeros(2*(n-1),1);

for i = 1:(n-1)
    
    A(i,permone(i)) = 1;
    A(i,permone(i+1)) = -1;
    
end

for i = 1:(n-1)
    
    A(i+(n-1),n + permtwo(i)) = 1;
    A(i+(n-1),n + permtwo(i+1)) = -1;
    
end

Aeq = zeros(2,2*n+1);
beq = zeros(2,1);


for i = 1:n
    
    Aeq (1,i) = 1;
    Aeq(2,i+n) = 1;
    
end



options = optimset('Algorithm','interior-point','MaxFunEvals', 30000,'MaxIter',10000);







[x,fval,exitflag] = fmincon(@nestedfun,start_vector,A,b,Aeq,beq,...
    [],[],[],options);
  

   
        for i = 1:n
            
            coordone(i) = x(i);
            coordtwo(i) = x(i+n);
        end
        
        addcon = x(2*n+1);
        
        meanprox = (1/(n*n))*sum(sum(prox));
        
        denom = 0;
        
        for i = 1:n
            for j = (i+1):n
                
                denom = denom+ (prox(i,j) - meanprox)^2;
                
            end
        end
        
        vaf = 1 - (fval/denom);
        
         function y = nestedfun(x)
        
        y = 0;
        
        for i = 1:n
            for j = (i+1):n
                
                y = y + (prox(i,j) - (x(2*n+1) + ((abs(x(i)-x(j)))^p + ...
                    (abs(x(i+n) - x(j+n)))^p)^(1/p)))^2;
                
            end
        end
         end
        

     fitted = zeros(n,n);
     
     for i = 1:n
         for j = 1:n
             
             if (i ~= j)
                 
                 fitted(i,j) = x(2*n+1) + ((abs(x(i)-x(j)))^p + ...
                    (abs(x(i+n) - x(j+n)))^p)^(1/p);
                
             end
             
         end
         
     end
                
        
   plot(coordone,coordtwo, 'k.')
   
   hold on
   
   for i = 1:n
       
       objectlabels{i,1} = int2str(i);
       
   end
   
   text(coordone,coordtwo,objectlabels,'fontsize', 10, ...
       'verticalalignment','bottom')
                
                




end

