############################ mpiBox.txt ##########################
#  toolbox mpiBox.txt for Eigenmath.           (2020) Dr.W.Lindner  
#  Use with the script    Moore-Penrose-pseudoinverse-Greville.pdf
############################ mpiBox.txt ##########################


mpi1(A)=do(At=transpose(A),         -- if A is matrix with detA<>0
           dot(inv(dot(At,A)),At) )

--------------------------- best fits ----------------------------

bestFit(B,A)= dot(mpi1(A),B)                             -- the Xo

proj1(B,A) =  dot(A, mpi1(A), B) 
                             -- orthogonal projection of B onto 
                             -- column space of A: the proximum Bo
            
oProj(B,A) = dot(A, dot( dot(A,A)^(-1), A, B)) 
             

EM(i,j,k,d) = do( M=unit(d,d), M[i,j]=k, M)

mpi(A) = test( rank(A)=1,              --if A vector & dot(A,A)<>0
               A * 1/dot(A,A) ,
               rank(A)>1,              --if A matrix & det(A)<>0
               do(At=transpose(A), dot(inv(dot(At,A)),At)) )
         
proj(B,A) = test( 
              rank(A)=1,               --if A vector & dot(A,A)<>0
              dot(A,B,A) / dot(A,A) ,
              rank(A)>1,               --if A matrix & det(A)<>0
              dot(A, mpi(A), B) )
               

Xfit(B,A) = dot(mpi(A),B)              -- the best fit solution Xo

Yfit(B,A) = dot(A, mpi(A), B)          -- the proximum Bo

---------------------------------------------------------
-- solution set of Linear System AX=B using pinv P
---------------------------------------------------------
isPinv(P,A)= test(dot(A,P,A) == A, 
                 "is pseudoinvers", 
                 "is NOT pseudoinvers")

isSolvable(A,B,P) = test(and(dot(A,P,A)==A, dot(A,P,B)==B), 
                        "is solvable", 
                        "is NOT solvable")

solSet(A,B,P,X) =  do(        -- solution set of LinearSystem AX=B
                      m=dim(dot(P,A),1), 
                      n=dim(dot(P,A),2),
                      E1=unit(m,n),
                      dot(P,B) + dot( (E1-dot(P,A)), X)     --LSSS
                     )

isMPI(P,A)= test(and(
                     dot(A,P,A)==A,
                     dot(P,A,P)==P,
                     transpose(dot(P,A))==dot(P,A),
                     transpose(dot(A,P))==dot(A,P) ),
                 " is MPI",
                 " is NOT the MPI")

mpi0(A)= dot( inv( dot(transpose(A),A) +0.0001*unit(dim(A,2)) ), transpose(A) ) 

mpiV(A) = test(  dot(A,A)==0,  
                 0*A,
                 1/dot(A,A) * A)

-----------    sequence of submatrices A1, A2, .. Ak i.e. A[:,1..k] 
Ai(k) =  test( k=1, a[1],                     -- case ai is vector
           do( AA=zero(k,dim(A,1)),           -- else k > 1 
               for(i,1,k, AA[i]=a[i]),        
               transpose(AA) ))

----------------- Greville m x n ------------------   A  typ m x n
--      procedure  GREVILLE  : A1+, A2+, .., An+  =   A+
---------------------------------------------------   A+ typ n x m

Greville(A)=
     do( n=dim(A,2) ,        
         m=dim(A,1) ,        
         a = transpose(A) ,
         null = zero(2,m) , 
       do(  -- k=1
         Ap = null ,  
         Ap[1] = mpiV(a[1]) ,
            -- k=2
         di = dot(Ap[1],a[2]) ,
         c  = a[2] - Ai(1)*di  ,
         test( c == null[1] ,
               b = (1+dot(di,di))^(-1)*dot(di,Ap[1]),
               b = mpiV(c)) ,
         B  = Ap[1] - di*b ,     
         Ap = zero(2,m)  ,
               Ap[1] = B ,
               Ap[2] = b,  
            -- k > 2
         do(for(k,3,n,
             di = dot(Ap,a[k]) ,   
             c  = a[k] - dot(Ai(k-1),di),
             test( c == null[1],
                   b = (1+dot(di,di))^(-1)*dot(di,Ap),
                   b = mpiV(c) ),
             B = Ap - outer(di,b),
             -- print(B),                    
             Ap = zero(k,m),
                  for(i,1,k-1, Ap[i] = B[i] ), 
                  Ap[k] = b, 
             Ap), -- close for, WATCH THE','
         Ap) -- close do around for-loop    
        )   -- close inner do
       )  -- close outer do               
    
-----------------  END GREVILLE -----------------   
############################ END mpiBox.txt ######################