function linear_regression_CLS
%
% linear_regression_CLS is a function that fits a straight line y = b*x + c
%   to a set of x,y data points using Combined Least Squares that allows 
%   for proper consideration of variances and covariances of the x,y data.
%   The function reads an ASCII textfile containing coordinate pairs (x,y)
%   standard deviations (sx,sy) and covariance (sxy).  In addition point 
%   name is included in the textfile.
%   Part of an ASCII textfile is shown below
%
% Data for Example 2
%
% Serial    x       sd_x           y      sd_y        covariance
%   1      -40.0    1.414213562   -24.0   1.732050808     0.5
%   2      -15.0    2.828427125   -24.0   2.236067978    -4.0
%   3       10.0    2.828427125   -12.0   2.645751311    -3.0
%   4       38.0    1.0            15.0   1.414213562     0.5
%   5       67.0    2.449489743    30.0   3.464101615     1.0 
%
%   Results are written to a textfile having the same path and name as the 
%   data file but with the extension ".out"

%============================================================================
% Function:  linear_regression_CLS
%
% Author:
%  Rod Deakin, 
%  BONBEACH VIC 3196
%  AUSTRALIA
%  email: randm.deakin@gmail.com
%
% Date:
%  Version 1.0  10 November 2014
%  Version 1.1  12 August   2015
%
% Remarks:
%  This function reads numeric data from a textfile containing coordinate
%  pairs (x,y) standard deviation and covariances associated with each pair 
%  and computes the parameters m and c of a line of best fit y = b*x + c 
%  using Combined Least Squares.  This allows proper consideration of the
%  standard deviations and covariances of the measurements (x,y coords).
%  Results are written to the screen and to textfile having the same path 
%  and name as the data file but with the extension ".out"
%
% Arrays:
%  A       - (m,n) coeff matrix of observation equations A*v + B*x = f
%  B       - (m,u) coeff matrix of observation equations A*v + B*x = f
%  f       - (m,1) vector of numeric terms
%  k       - (m,1) vector of Lagrange multipliers k = We*(f-B*x)
%  N       - (u,u) coefficient matrix of Normal equations N*x = t and 
%                  N = Bt*We*B
%  Ninv    - (u,u) inverse of coefficient matrix of Normal equations
%  Q       - (n,n) matrix of estimates of variances and covariances
%  Qll     - (n,n) cofactor matrix of adjusted observations
%  Qvv     - (n,n) cofactor matrix of residuals
%  We      - (m,m) equivalent weight matrix We = inv(A*Q*At)
%  x_coord - (m,1) vector of x coordinates 
%  y_coord - (m,1) vector of y coordinates
%  serial  - (m,1) vector of station serial numbers
%  sx      - (m,1) vector of st. dev's of x coordinates
%  sxy     - (m,1) vector of covariances between x and y coordinates
%  sy      - (m,1) vector of st. dev's of y coordinates
%  t       - (u,1) vector of numeric terms of Normal equations N*x = t and
%                  t = Bt*We*f
%  v       - (n,1) vector of residuals v = Q*At*k
%  W       - (n,n) matrix of weight coefficients W = inv(Q)
%  x       - (u,1) vector of solutions x = Ninv*t
%
% Variables
%  b       - b = tan(theta) = (y2-y1)/(x2-x1) = slope of straight line
%  b0      - approximate value of slope of straight line
%  c       - point on y-axis where straight line intersects
%  db      - small correction to slope of line b = b0 + db
%  eps     - small quantity for testing for convergence
%  iter    - iteration count
%  m       - number of observation equations equal to number of (x,y) pairs
%  n       - number of observations (or measurements)
%  u       - number of unknowns (parameters m and c)
%  i,j,kk  - integer counters
%  vfact   - variance factor
%  xc,yc   - x and y centroids (averages)
%  x1,x2   - x coords of first and last data points
%  y1,y2   - y coords of first and last data points
%
% References:
%  Deakin, R.E., 2014, 'Linear Regression Using Combined Least Squares',
%     Private Notes, Bonbeach VIC 3196, Version 1, 03-Nov-2014
%
%============================================================================

%-------------------------------------------------------------------------
% 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);
serial  = fileTemp(:,1);
x_coord = fileTemp(:,2);
sx      = fileTemp(:,3);
y_coord = fileTemp(:,4);
sy      = fileTemp(:,5);
sxy     = fileTemp(:,6);
%--------------------------------------------------------------------------
% Determine the adjustment constants m, n and u
%--------------------------------------------------------------------------
% determine the number of observation equations 
m = length(x_coord);
% set the number of measurements
n = 2*m;
% set the number of unknowns
u = 2;
%--------------------------------------------------------------------------
% Print heading and given data
%--------------------------------------------------------------------------
% open the output file print the data
fidout  = fopen(outfilepath,'wt');
fprintf(fidout,'\n\n Linear Regression Using Combined Least Squares');
fprintf(fidout,'\n ==============================================');
fprintf(fidout,'\n solution for the parameters b and c in the equation for the line of best fit');
fprintf(fidout,'\n y = b*x + c  where b = tan(theta) is the slope and c is the intercept on the y-axis');
fprintf(fidout,'\n\n Input Data');
fprintf(fidout,'\n Serial      x       sx              y       sy            sxy');
for kk = 1:m
    fprintf(fidout,'\n %4d %12.4f  %8.6f  %12.4f  %8.6f     %9.6f',...
            serial(kk,1),x_coord(kk,1),sx(kk,1),y_coord(kk,1),...
            sy(kk,1),sxy(kk,1));
end
fprintf(fidout,'\n\n there are n = %2d measurements (the x,y coordinates),',n); 
fprintf(fidout,'\n m = %2d observation equations and u = %2d unknowns in the',m,u); 
fprintf(fidout,'\n system of equations A*v + B*x = f'); 
%--------------------------------------------------------------------------
% set the elements of the cofactor matrix of observations Q
%--------------------------------------------------------------------------
Q = zeros(n,n);
for kk = 1:m
    j = 2*kk;
    Q(j-1,j-1) = sx(kk,1)^2;  % variance x = (sx)^2
    Q(j-1,j) = sxy(kk,1);     % covariance
    Q(j,j-1) = Q(j-1,j);      % covariance
    Q(j,j) = sy(kk,1)^2;      % variance y = (sy)^2
end
%--------------------------------------------------------------------------
% Reduce the coordinates to a centroidal system
%--------------------------------------------------------------------------
xc = sum(x_coord)/m;
yc = sum(y_coord)/m;
x_coord = x_coord-xc;
y_coord = y_coord-yc;
%--------------------------------------------------------------------------
% form the coefficient matrix B of the observation equation A*v +B*x = f
%--------------------------------------------------------------------------
B = zeros(m,u);
for kk = 1:m
    B(kk,1) = -x_coord(kk,1);
    B(kk,2) = -1.0;
end
%--------------------------------------------------------------------------
% set approximate value for slope b0 where b = b0 + db and db is a small
% correction.  b = tan(theta) = (y2-y1)/(x2-x1)
%--------------------------------------------------------------------------
y1 = y_coord(1,1);
y2 = y_coord(m,1);
x1 = x_coord(1,1);
x2 = x_coord(m,1);
b0 = (y2-y1)/(x2-x1);

%--------------------------------------------------------------------------
% Solution of least squares estimates by iteration
%--------------------------------------------------------------------------
% while loop for the computation of solution vector x where x(1,1) = db
% and x(2,1) = c.  Break from while loop when abs(db) < eps 
% set the smmal increment to the slope where b = b0 + db
db = 10;
% set the iteration number
iter = 1;
% set eps, a small quantity
eps = 1.0e-8;
while 1
    iter;
    if abs(db)<eps
        break
    end
    if iter > 10
        break
    end
    % form the coefficient matrix A of observation equations A*v + B*x = f
    A = zeros(m,n);
    for kk = 1:m
        j = 2*kk;
        A(kk,j)   = 1.0;
        A(kk,j-1) = -b0;
    end
    % form the vector of numeric terms f
    f = zeros(m,1);
    for kk = 1:m
        f(kk,1) = b0*x_coord(kk,1)-y_coord(kk,1);
    end
    % form ther equivalent weight matrix We
    We = inv(A*Q*A');
    % form the normal equation coefficient matrix N
    N = B'*We*B;
    % form the vector of numeric terms t
    t = B'*We*f;
    % solve the system Nx = t for the unknown parameters x
    Ninv = inv(N);
    x = Ninv*t;
    
    if iter == 1
        % print matrices A, B, f, Q, We, N, t and x
        fprintf(fidout,'\n\n System of equations and solution for first iteration');
        fprintf(fidout,'\n where approximate slope b0 = % 12.9f',b0);
        fprintf(fidout,'\n Coordinates have been reduced to a centroidal system');
        fprintf(fidout,'\n The * symbol indicates new row of a matrix');
        fprintf(fidout,'\n Coefficient matrix A');
        for i = 1:m
            kk = 1;
            fprintf(fidout,'\n* ');
            for j = 1:n
                if (kk > 6)
                    fprintf(fidout,'\n  ');
                    kk = 1;
                end
                fprintf(fidout,'% 12.9f ',A(i,j));
                kk = kk+1;
            end
        end
        fprintf(fidout,'\n Coefficient matrix B');
        for i = 1:m
            kk = 1;
            fprintf(fidout,'\n* ');
            for j = 1:u
                if (kk > 6)
                    fprintf(fidout,'\n  ');
                    kk = 1;
                end
                fprintf(fidout,'% 12.9f ',B(i,j));
                kk = kk+1;
            end
        end
        fprintf(fidout,'\n Vector f-transposed');
        kk = 1;
        fprintf(fidout,'\n  ');
        for i = 1:m
            if (kk > 6)
                fprintf(fidout,'\n  ');
                kk = 1;
            end
            fprintf(fidout,'% 12.9f ',f(i,1));
            kk = kk+1;
        end
        fprintf(fidout,'\n Cofactor matrix Q');
        for i = 1:n
            kk = 1;
            fprintf(fidout,'\n* ');
            for j = 1:n
                if (kk > 6)
                    fprintf(fidout,'\n  ');
                    kk = 1;
                end
                fprintf(fidout,'% 12.9f ',Q(i,j));
                kk = kk+1;
            end
        end
        fprintf(fidout,'\n Equivalent weight matrix We = inv(A*Q*At)');
        for i = 1:m
            kk = 1;
            fprintf(fidout,'\n* ');
            for j = 1:m
                if (kk > 6)
                    fprintf(fidout,'\n  ');
                    kk = 1;
                end
                fprintf(fidout,'% 15.9f ',We(i,j));
                kk = kk+1;
            end
        end
        fprintf(fidout,'\n Matrix N = Bt*We*B');
        for i = 1:u
            kk = 1;
            fprintf(fidout,'\n* ');
            for j = 1:u
                if (kk > 6)
                    fprintf(fidout,'\n  ');
                    kk = 1;
                end
                fprintf(fidout,'% 18.9f ',N(i,j));
                kk = kk+1;
            end
        end
        fprintf(fidout,'\n Vector t-transposed where t = Bt*We*f');
        kk = 1;
        fprintf(fidout,'\n  ');
        for i = 1:u
            if (kk > 6)
                fprintf(fidout,'\n  ');
                kk = 1;
            end
            fprintf(fidout,'% 12.9f ',t(i,1));
            kk = kk+1;
        end
        fprintf(fidout,'\n Vector of solutions x-transposed where x = Ninv*t');
        kk = 1;
        fprintf(fidout,'\n  ');
        for i = 1:u
            if (kk > 6)
                fprintf(fidout,'\n  ');
                kk = 1;
            end
            fprintf(fidout,'% 12.9f ',x(i,1));
            kk = kk+1;
        end
    end
    % update value of b0
    db = x(1,1);
    b0 = b0+db;
    % increment iteration count
    iter = iter+1;
end
%--------------------------------------------------------------------------
% Set the solutions for slope b and y-axis intercept c for the straight
% line y = bx + c
%--------------------------------------------------------------------------
db = x(1,1);
b  = b0+db;
c  = b*(-xc)+x(2,1)+yc;
%--------------------------------------------------------------------------
% compute vector of Lagrange multipliers k and vector of residuals v
%--------------------------------------------------------------------------
k = We*(f-B*x);
v = Q*A'*k;

%--------------------------------------------------------------------------
% compute the variance factor and cofactor matrices Qxx, Qll, Qvv
%--------------------------------------------------------------------------
W = inv(Q);
vfact = (v'*W*v)/(m-u);

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

% compute the cofactor matrix of the adjusted quantities
Qll = Q + Q*A'*We*B*Ninv*B'*We*A*Q - Q*A'*We*A*Q;

% compute the cofactor matrix of the residuals
Qvv = Q - Qll;

%--------------------------------------------------------------------------
% Print solutions
%--------------------------------------------------------------------------
fprintf(fidout,'\n\n Solutions after %2d iterations',iter);
fprintf(fidout,'\n ==============================');
fprintf(fidout,'\n\n slope of line m = tan(theta) = % 8.6f',b');
fprintf(fidout,'\n y-axis intercept           c = % 8.6f',c');
fprintf(fidout,'\n\n Vector of Lagrange multipliers k-transposed where k = We*(f-B*x)');
kk = 1;
fprintf(fidout,'\n  ');
for i = 1:m
    if (kk > 6)
        fprintf(fidout,'\n  ');
        kk = 1;
    end
    fprintf(fidout,'% 12.9f ',k(i,1));
    kk = kk+1;
end
fprintf(fidout,'\n\n Vector of residuals v-transposed where v = Q*At*k');
kk = 1;
fprintf(fidout,'\n  ');
for i = 1:n
    if (kk > 6)
        fprintf(fidout,'\n  ');
        kk = 1;
    end
    fprintf(fidout,'% 12.9f ',v(i,1));
    kk = kk+1;
end

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

fprintf(fidout,'\n\n Cofactor matrix of unknowns Qxx (upper triangular part)');
fprintf(fidout,'\n The * symbol indicates new row beginning at diagonal');
for j = 1:u
  fprintf(fidout,'\n* ');
  count = 0;
  for kk = j:u
    if count > 5
      count = 0;
      fprintf(fidout,'\n  ');
    end  
    fprintf(fidout,'%14.6e',Qxx(j,kk));
    count = count + 1;
  end  
end  

fprintf(fidout,'\n\n Cofactor matrix of residuals Qvv (upper triangular part)');
fprintf(fidout,'\n The * symbol indicates new row beginning at diagonal');
for j = 1:n
  fprintf(fidout,'\n* ');
  count = 0;
  for kk = j:n
    if count > 5
      count = 0;
      fprintf(fidout,'\n  ');
    end  
    fprintf(fidout,'%14.6e',Qvv(j,kk));
    count = count + 1;
  end  
end  

fprintf(fidout,'\n\n Cofactor matrix of adjusted observations Qll (upper triangular part)');
fprintf(fidout,'\n The * symbol indicates new row beginning at diagonal');
for j = 1:n
  fprintf(fidout,'\n* ');
  count = 0;
  for kk = j:n
    if count > 5
      count = 0;
      fprintf(fidout,'\n  ');
    end  
    fprintf(fidout,'%14.6e',Qll(j,kk));
    count = count + 1;
  end  
end  

fprintf(fidout,'\n\n');
 
% close the output file
fclose(fidout);
