#  MAPLE.  October 2, 2017.
#
#  Implement the affine scaling method.
#
#  Input: A  - an m x n matrix with full row rank.
#         b  - an m x 1 matrix 
#         c  - and n x 1 matrix
#         x0 - a strictly positive solution to Ax=b
#         r  - a real number  0 < r < 1, the step size.
#
#  procedure "affineit" performs one iteration of the affine scaling method.
#  Output: x1 - an n x 1 matrix which is the next strictly positive 
#               feasible solution after x0 in the affine scaling method.
#  The code prints most intermediate calculations to stdout
#  Be careul when using sqrt(.), etc. since MAPLE cannot always tell
#  if entries are non-zero. It is simplest to use floating point values.

with(LinearAlgebra):



######################################
# Input arguments: matrix A, RHS b, objective c, current soln x0, step size r
######################################
affineit := proc( oA, b, c, x0, r)
local  A,cc,  n,m,x,X0,i,tA,P,d,beta,j,ones;
# Assume: Given A, b, c with A (m by n) having full row rank and 
#   a strictly positive solution x0 (n by 1).
n := ColumnDimension(oA);
m := RowDimension(oA);

x := x0;


   print("---------------------------------------------------");
   print("Beginning of iteration");
   print("---------------------------------------------------");

   # Build X0
   X0 := Matrix(n,n,0):
   for i to n do X0[i,i] := x[i,1]; od:
   print("   X0 = ",X0 );
   # Build A-hat and c-hat
   A := oA.X0;  # oA is the "original" matrix A
   cc :=X0.c;
   print("   A-hat = ",A ,"   c-hat = ",cc );
   # Compute projection matrix P
   tA := Transpose(A);
   P := IdentityMatrix(n) - tA.( A . tA )^(-1) . A;
   print(" A * A^T = ", A . tA );
   print(" (A * A^T)^(-1) = ", MatrixInverse( A.tA)  );
   print(" A^T * (A * A^T)^(-1) = ", tA.MatrixInverse( A . tA)  );
   print(" A^T * (A * A^T)^(-1) * A = ",  tA.MatrixInverse( A . tA).A  );
   print(" P = ", P );
   # Find search direction
   d := P.cc;
   print( " d = ", d );
   # Compute ratios
   beta := infinity;
   for j to n do if d[j,1] < 0 then beta := min( beta, -1/d[j,1] ); fi od:
   print(" beta = ",beta);
   # Scaled update
   ones := Matrix(n,1,1.0);  # all ones vector
   x := ones + r*beta*d ;
   print(" new scaled solution x-hat = X0 x = ",x );
   # Next iterate is obtained by scaling back to original problem.
   x := X0.x;
   print(" new solution  x = ", x ," objective value zeta =", 
      Transpose(c).x );

RETURN( x );
end;