function least_squares
%
% LEAST_SQUARES reads an ASCII textfile containing data for a least squares
%  adjustment of indirect observations in the matrix form v + Bx = f where
%  v is an n,1 vector of residuals, B is an n,u coefficient matrix, x is a
%  u,1 vector of unknowns and f is an n,1 vector of constant terms.  n is
%  the number of equations and u is the number of unknowns.  The data textfile
%  contains n rows of data with u+2 columns.  The elements of B are
%  contained in the first u columns.  The numeric terms (f) is contained in
%  column u+1 and the weight of each observation is contained in column
%  u+2.  Results are written to a textfile having the same path and name as 
%  the data file but with the extension ".out" 

%============================================================================
% Function:  least_squares
%
% Author:
%  Rod Deakin, 
%  Department of Geospatial Science, RMIT University,
%  GPO Box 2476V, MELBOURNE VIC 3001
%  AUSTRALIA
%  email: randm.deakin@gmail.com
%
% Date:
%  Version 1.0  14 April  2003
%  Version 1.1   9 June   2003
%  Version 1.1a 12 August 2015 (Change of email address)
%
% Remarks:
%  This function reads numeric data from a textfile containing data for a
%  least squares adjustment of indirect observations in the form of 
%  v + Bx = f where v is an n,1 vector of residuals, B is an n,u matrix of
%  coefficients, x in a u,1 vector of unknown parameters and f is an n,1
%  vector of numeric terms.  n is the number of equations and u is the 
%  number of unknowns.  The data text file contains n rows of data in u+2
%  columns.  Column u+1 are the numeric terms, column u+2 are the weights
%  of each observation.  Results are written to a textfile having the same
%  path name as the data file but with the extension ".out"
%
% Arrays:
%  B       - (n,u) coeff matrix of observation equation v + Bx = f
%  f       - (n,1) vector of numeric terms
%  N       - (u,u) coefficient matrix of Normal equations Nx = t
%  Ninv    - (u,u) inverse of N
%  Qxx     - (u,u) cofactor matrix of unknowns
%  Qvv     - (n,n) cofactor matrix of residuals
%  Qll     - (n,n) cofactor matrix of adjusted measurements
%  t       - (u,1) vector of numeric terms of Normal equations Nx = t
%  v       - (n,1) vector of residuals
%  W       - (n,n) weight matrix
%  w       - (n,1) vector of weights
%  x       - (u,1) vector of solutions
%
% Variables
%  n       - number of equations
%  u       - number of unknowns
%  vfact   - variance factor
%
% References:
%  Notes on Least Squares (2003), Department of Geospatial Science, RMIT
%  University, 2003
%
%============================================================================

%-------------------------------------------------------------------------
% 1. Call the User Interface (UI) to choose the input data file name
% 2. Concatenate strings to give the path and file name of the input file
% 3. Strip off the extension from the file name to give the rootName
% 4. Add extension ".out" to rootName to give the output filename
% 5. Concatenate strings to give the path and file name of the output file
%-------------------------------------------------------------------------
filepath = strcat('d:\temp\','*.txt');
[infilename,inpathname] = uigetfile(filepath);
infilepath = strcat(inpathname,infilename);
rootName   = strtok(infilename,'.');
outfilename = strcat(rootName,'.out');
outfilepath = strcat(inpathname,outfilename);

%----------------------------------------------------------
% 1. Load the data into an array whose name is the rootName
% 2. set fileTemp = rootName
% 3. Copy columns of data into individual arrays
%----------------------------------------------------------
load(infilepath);
fileTemp = eval(rootName);
% get the number of rows (n) and the number of columns (m)
[n,m] = size(fileTemp);
% set the number of unknowns
u = m-2;
% copy the data into B, f and w
B = fileTemp(:,1:u);
f = fileTemp(:,u+1);
w = fileTemp(:,m);

% set the elements of the weight matrix W
W = zeros(n,n);
for k = 1:n
  W(k,k) = w(k);
end  

% form the normal equation coefficient matrix N
N = B'*W*B;

% form the vector of numeric terms t
t = B'*W*f;

% solve the system Nx = t for the unknown parameters x
Ninv = inv(N);
x = Ninv*t;

% compute residuals
v = f - (B*x);

% compute the variance factor
vfact = (f'*W*f - x'*t)/(n-u);

% compute the cofactor matrix of the adjusted quantities
Qxx = Ninv;

% compute the cofactor matrix of the residuals
Qvv = inv(W)-B*Ninv*B';

% compute the cofactor matrix of the adjusted quantities
Qll = inv(W)-Qvv;

% open the output file print the data
fidout  = fopen(outfilepath,'wt');

fprintf(fidout,'\n\nLeast Squares Adjustment of Indirect Observations');
fprintf(fidout,'\n\nInput Data');
fprintf(fidout,'\n\nCoefficient matrix B of observation equations v + Bx = f');
for j = 1:n
  fprintf(fidout,'\n');
  for k = 1:u
    fprintf(fidout,'%15.6f',B(j,k));
  end  
end  

fprintf(fidout,'\n\nVector of numeric terms f and weights w of observation equations v + Bx = f');
for k = 1:n
  fprintf(fidout,'\n%15.6f  %15.6f',f(k,1),w(k));
end  

fprintf(fidout,'\n\nCoefficient matrix N of Normal equations Nx = t');
fprintf(fidout,'\n(upper triangular part)');
for j = 1:u
  fprintf(fidout,'\n');
  count = 0;
  for k = j:u
    if count > 4
      count = 0;
      fprintf(fidout,'\n');
    end  
    fprintf(fidout,'%15.6f',N(j,k));
    count = count + 1;
  end  
end  

fprintf(fidout,'\n\nVector of numeric terms t of Normal equations Nx = t');
for k = 1:u
  fprintf(fidout,'\n%15.6f',t(k,1));
end  

fprintf(fidout,'\n\nInverse of Normal equation coefficient matrix');
fprintf(fidout,'\n(upper triangular part)');
for j = 1:u
  fprintf(fidout,'\n');
  count = 0;
  for k = j:u
    if count > 4
      count = 0;
      fprintf(fidout,'\n');
    end  
    fprintf(fidout,'%16.4e',Ninv(j,k));
    count = count + 1;
  end  
end  

fprintf(fidout,'\n\nVector of solutions x');
for k = 1:u
  fprintf(fidout,'\n%20.9e',x(k,1));
end  

fprintf(fidout,'\n\nVector of residuals v');
for k = 1:n
  fprintf(fidout,'\n%12.6f',v(k,1));
end

fprintf(fidout,'\n\nVariance factor = %10.6e',vfact);

fprintf(fidout,'\n\nCofactor matrix of unknowns Qxx');
fprintf(fidout,'\n(upper triangular part)');
for j = 1:u
  fprintf(fidout,'\n');
  count = 0;
  for k = j:u
    if count > 4
      count = 0;
      fprintf(fidout,'\n');
    end  
    fprintf(fidout,'%16.4e',Qxx(j,k));
    count = count + 1;
  end  
end  

fprintf(fidout,'\n\nCofactor matrix of residuals Qvv');
fprintf(fidout,'\n(upper triangular part)');
for j = 1:n
  fprintf(fidout,'\n');
  count = 0;
  for k = j:n
    if count > 4
      count = 0;
      fprintf(fidout,'\n');
    end  
    fprintf(fidout,'%16.4e',Qvv(j,k));
    count = count + 1;
  end  
end  

fprintf(fidout,'\n\nCofactor matrix of adjusted observations Qll');
fprintf(fidout,'\n(upper triangular part)');
for j = 1:n
  fprintf(fidout,'\n');
  count = 0;
  for k = j:n
    if count > 4
      count = 0;
      fprintf(fidout,'\n');
    end  
    fprintf(fidout,'%16.4e',Qll(j,k));
    count = count + 1;
  end  
end  

fprintf(fidout,'\n\n');

% close the output file
fclose(fidout);
