' ******************************************************************
' Correlation and principal component analysis
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE MatOk 0
' No error

#DEFINE MatNonConv -1
' Non-convergence

#DEFINE MatErrDim -3
' Non-compatible dimensions

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the last function evaluation

' ------------------------------------------------------------------
' External function
' ------------------------------------------------------------------

DECLARE SUB EigenSym(A()      AS DOUBLE,_
                     Lambda() AS DOUBLE,_
                     V()      AS DOUBLE)
' Eigenvalues / eigenvectors of a positive semi-definite symmetric matrix

' ******************************************************************

SUB VecMean(X() AS DOUBLE, M() AS DOUBLE)
' ------------------------------------------------------------------
' Computes the mean vector M() from matrix X()
' ------------------------------------------------------------------

  DIM AS INTEGER Lb1, Ub1, Lb2, Ub2, I, J, N
  DIM AS DOUBLE  Sum

  Lb1 = LBOUND(X, 1)
  Ub1 = UBOUND(X, 1)
  Lb2 = LBOUND(X, 2)
  Ub2 = UBOUND(X, 2)

  IF LBOUND(M) <> Lb2 OR UBOUND(M) <> Ub2 THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  ErrCode = MatOk
  N = Ub1 - Lb1 + 1

  FOR J = Lb2 TO Ub2
    Sum = 0
    FOR I = Lb1 TO Ub1
      Sum = Sum + X(I, J)
    NEXT I
    M(J) = Sum / N
  NEXT J
END SUB

SUB VecSD(X() AS DOUBLE, M() AS DOUBLE, S() AS DOUBLE)
' ------------------------------------------------------------------
' Computes the vector of standard deviations S() from matrix X()
' and mean vector M()
' ------------------------------------------------------------------

  DIM AS INTEGER Lb1, Ub1, Lb2, Ub2, I, J, N
  DIM AS DOUBLE  Sum, D

  Lb1 = LBOUND(X, 1)
  Ub1 = UBOUND(X, 1)
  Lb2 = LBOUND(X, 2)
  Ub2 = UBOUND(X, 2)

  IF LBOUND(M) <> Lb2 OR UBOUND(M) <> Ub2 OR _
     LBOUND(S) <> Lb2 OR UBOUND(S) <> Ub2    _
  THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  ErrCode = MatOk
  N = Ub1 - Lb1 + 1

  FOR J = Lb2 TO Ub2
    Sum = 0
    FOR I = Lb1 TO Ub1
      D = X(I,J) - M(J)
      Sum = Sum + D * D
    NEXT I
    S(J) = SQR(Sum / N)
  NEXT J
END SUB

SUB MatVarCov(X() AS DOUBLE, M() AS DOUBLE, V() AS DOUBLE)
' ------------------------------------------------------------------
' Computes the variance-covariance matrix V() from matrix X()
' and mean vector M()
' ------------------------------------------------------------------

  DIM AS INTEGER Lb1, Ub1, Lb2, Ub2, I, J, K, N
  DIM AS DOUBLE  Sum

  Lb1 = LBOUND(X, 1)
  Ub1 = UBOUND(X, 1)
  Lb2 = LBOUND(X, 2)
  Ub2 = UBOUND(X, 2)

  IF LBOUND(M)    <> Lb2 OR UBOUND(M)    <> Ub2 OR _
     LBOUND(V, 1) <> Lb2 OR UBOUND(V, 1) <> Ub2 OR _
     LBOUND(V, 2) <> Lb2 OR UBOUND(V, 2) <> Ub2    _
  THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  ErrCode = MatOk
  N = Ub1 - Lb1 + 1

  FOR I = Lb2 TO Ub2
    FOR J = I TO Ub2
      Sum = 0
      FOR K = Lb1 TO Ub1
        Sum = Sum + (X(K,I) - M(I)) * (X(K,J) - M(J))
      NEXT K
      V(I,J) = Sum / N
      IF I <> J THEN V(J,I) = V(I,J)
    NEXT J
  NEXT I
END SUB

SUB MatCorrel(V() AS DOUBLE, R() AS DOUBLE)
' ------------------------------------------------------------------
' Computes the correlation matrix R() from the variance-covariance
' matrix V()
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, I, J

  Lb = LBOUND(V, 1)
  Ub = UBOUND(V, 1)

  IF LBOUND(V, 2) <> Lb OR UBOUND(V, 2) <> Ub OR _
     LBOUND(R, 1) <> Lb OR UBOUND(R, 1) <> Ub OR _
     LBOUND(R, 2) <> Lb OR UBOUND(R, 2) <> Ub    _
  THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  ErrCode = MatOk

  FOR I = Lb TO Ub
    R(I,I) = 1
    FOR J = I + 1 TO Ub
      R(I,J) = V(I,J) / SQR(V(I,I) * V(J,J))
      R(J,I) = R(I,J)
    NEXT J
  NEXT I
END SUB

SUB PCA(R() AS DOUBLE, Lambda() AS DOUBLE, C() AS DOUBLE, Rc() AS DOUBLE)
' ------------------------------------------------------------------
' Performs a principal component analysis of a correlation matrix
' ------------------------------------------------------------------
' Input : R()      = Correlation matrix
' Output: Lambda() = Eigenvalues of the correlation matrix
'                    (in descending order)
'         C()      = Eigenvectors of the correlation matrix
'                    (column J contains the J-th eigenvector)
'         Rc()     = Correlations between principal factors and
'                    variables. Rc(I,J) is the correlation
'                    coefficient between variable I and factor J.
' ------------------------------------------------------------------
' Possible results : MatOk      : No error
'                    MatNonConv : Non-convergence of eigenvalue
'                                    determination
' ------------------------------------------------------------------
' NB : This sub destroys the original matrix R
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, I, J
  DIM AS DOUBLE  Rac

  Lb = LBOUND(Lambda)
  Ub = UBOUND(Lambda)

  IF LBOUND(R, 1) <> Lb OR UBOUND(R, 1) <> Ub OR _
     LBOUND(R, 2) <> Lb OR UBOUND(R, 2) <> Ub OR _
     LBOUND(C, 1) <> Lb OR UBOUND(C, 1) <> Ub OR _
     LBOUND(C, 2) <> Lb OR UBOUND(C, 2) <> Ub OR _
     LBOUND(Rc,1) <> Lb OR UBOUND(Rc,1) <> Ub OR _
     LBOUND(Rc,2) <> Lb OR UBOUND(Rc,2) <> Ub    _
  THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  ' Compute eigenvalues and eigenvectors of correlation matrix
  EigenSym R(), Lambda(), C()

  IF ErrCode <> MatOk THEN EXIT SUB

  ' Compute correlations between principal factors and reduced variables
  FOR J = Lb TO Ub
    Rac = SQR(Lambda(J))
    FOR I = Lb TO Ub
      Rc(I,J) = C(I,J) * Rac
    NEXT I
  NEXT J
END SUB

SUB ScaleVar(X() AS DOUBLE, M() AS DOUBLE, S() AS DOUBLE, Z() AS DOUBLE)
' ------------------------------------------------------------------
' Scales a set of variables by subtracting means and dividing by SD
' ------------------------------------------------------------------
' Input : X() = matrix original variables (columns)
'         M() = mean vector
'         S() = vector of standard deviations
' Output: Z() = matrix of scaled variables
' ------------------------------------------------------------------

  DIM AS INTEGER Lb1, Ub1, Lb2, Ub2, I, J

  Lb1 = LBOUND(X, 1)
  Ub1 = UBOUND(X, 1)
  Lb2 = LBOUND(X, 2)
  Ub2 = UBOUND(X, 2)

  IF LBOUND(M)    <> Lb2 OR UBOUND(M)    <> Ub2 OR _
     LBOUND(S)    <> Lb2 OR UBOUND(S)    <> Ub2 OR _
     LBOUND(Z, 1) <> Lb1 OR UBOUND(Z, 1) <> Ub1 OR _
     LBOUND(Z, 2) <> Lb2 OR UBOUND(Z, 2) <> Ub2    _
  THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  ErrCode = MatOk

  FOR I = Lb1 TO Ub1
    FOR J = Lb2 TO Ub2
      Z(I,J) = (X(I,J) - M(J)) / S(J)
    NEXT J
  NEXT I
END SUB

SUB PrinFac(Z() AS DOUBLE, C() AS DOUBLE, F() AS DOUBLE)
' ------------------------------------------------------------------
' Computes principal factors
' ------------------------------------------------------------------
' Input : Z() = matrix of scaled variables (columns)
'         C() = matrix of eigenvectors from PCA
' Output: F() = matrix of principal factors
' ------------------------------------------------------------------

  DIM AS INTEGER Lb1, Ub1, Lb2, Ub2, I, J, K

  Lb1 = LBOUND(Z, 1)
  Ub1 = UBOUND(Z, 1)
  Lb2 = LBOUND(Z, 2)
  Ub2 = UBOUND(Z, 2)

  IF LBOUND(C, 1) <> Lb2 OR UBOUND(C, 1) <> Ub2 OR _
     LBOUND(C, 2) <> Lb2 OR UBOUND(C, 2) <> Ub2 OR _
     LBOUND(F, 1) <> Lb1 OR UBOUND(F, 1) <> Ub1 OR _
     LBOUND(F, 2) <> Lb2 OR UBOUND(F, 2) <> Ub2    _
  THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  ErrCode = MatOk

  FOR I = Lb1 TO Ub1
    FOR J = Lb2 TO Ub2
      F(I,J) = 0
      FOR K = Lb2 TO Ub2
        F(I,J) = F(I,J) + Z(I,K) * C(K,J)
      NEXT K
    NEXT J
  NEXT I
END SUB

