Matrix and Statistical operations

lib/plugins/matrix/matrix.c

@@ -419,7 +419,7 @@ PROCEDURE(mprint) {
         printf("Matrix validation for %d failed: %d\n",GETINT(ARG0),status);
         RETURNINT(status);
     }
-    printMatrix(GETINT(ARG0),"");
+    printMatrix(GETINT(ARG0),GETSTRING(ARG1));
     ENDPROC
 }
 
@@ -2428,7 +2428,7 @@ LOADFUNCS
     ADDPROC(mset,      "matrix.mset",      "b",  ".int", "m0=.int, row=.int, col=.int,value=.float");
     ADDPROC(mget,      "matrix.mget",      "b",  ".float","m0=.int, row=.int, col=.int");
     ADDPROC(mexpand,   "matrix.mexpand",   "b",  ".int", "m0=.int, newrows=.int, init=.float");
-    ADDPROC(mprint,    "matrix.mprint",    "b",  ".int", "m0=.int");
+    ADDPROC(mprint,    "matrix.mprint",    "b",  ".int", "m0=.int,hdr=.string");
     ADDPROC(mfree,     "matrix.mfree",     "b",  ".int", "m0=.int");
     ADDPROC(mdet,      "matrix.mdet",      "b",  ".int", "m0=.int");
     ADDPROC(mlu,       "matrix.mlu",       "b",  ".int", "m0=.int, L=.string, U=.string");

lib/plugins/matrix/matrix.md

@@ -20,11 +20,8 @@ This library provides:
 - **Error Handling**: Comprehensive mechanisms to handle common matrix errors, ensuring robust and fail-safe operations.
 - **Performance Optimization**: Features for parallel processing, memory optimization, and numerical stability.
 
-The library is structured with best practices for memory management, error recovery, and performance tuning, making it suitable for both small-scale and high-performance computing environments. With its extensibility for advanced features like sparse matrices and system limits, the **Matrix Operations Library** delivers a powerful and flexible solution for handling complex matrix workflows.
-
 ---
 
-
 ## Basic Matrix Functions
 
 ### Matrix Creation
@@ -46,23 +43,23 @@ id = mcreate(3, 3, "myMatrix");
     - `m0`: First input matrix (int)
     - `m1`: Second input matrix (int)
     - `mid`: Result matrix identifier (string)
-- **Returns**: Matrix handle (int)
+- **Returns**: new Matrix handle (int)
 
 ### Scalar Multiplication (`mprod`,'description')
 - **Purpose**: Multiplies a matrix by a scalar: \( B = A * scalar \)
 - **Parameters**:
     - `m0`: Input matrix (int)
     - `prod`: Scalar value (float)
     - `mid`: Result matrix identifier (string)
-- **Returns**: Matrix handle (int)
+- **Returns**: new Matrix handle (int)
 
 
 ### Matrix Inversion (`minvert`,'description')
 - **Purpose**: Computes the inverse of a matrix: \( B = A^{-1} \)
 - **Parameters**:
     - `m0`: Input matrix (int)
     - `mid`: Result matrix identifier (string)
-- **Returns**: Matrix handle (int)
+- **Returns**: new Matrix handle (int)
 
 ### Matrix Transposition (`mtranspose`,'description')
 - **Purpose**: Computes the transpose of a matrix: \( B = A^T \)
@@ -82,22 +79,22 @@ mtrans(m1, "result"); // Transpose
 ### `mdet`
 - **Purpose**: Calculate matrix determinant
 - **Parameters**:
-  - m0: Input matrix (.int)
-- **Returns**: Success status (.int)
+  - m0: Input matrix (matrix-handle)
+- **Returns**: Success status (integer)
 
 ### `mrank`
 - **Purpose**: Calculate matrix rank
 - **Parameters**:
-  - m0: Input matrix (.int)
-- **Returns**: Matrix rank (.int)
+  - m0: Input matrix (matrix-handle)
+- **Returns**: Matrix rank (integer)
 ### `mlu`
 - **LU Decomposition**: `mlu(m0, L, U)` — Performs LU decomposition, returning lower (`L`) and upper (`U`) triangular matrices.
-- m0: Input matrix (.int)
+- m0: Input matrix (matrix-handle)
 ### `mprint`
 - **Purpose**: Print matrix to output
 - **Parameters**:
-  - m0: Matrix to print (.int,"alternative-heading") - if empty Matrix heading defined with MCREATE is used
-- **Returns**: Success status (.int)
+  - m0: Matrix to print (matrix-handle,"alternative-heading") - if empty Matrix heading defined with MCREATE is used
+- **Returns**: Success status (integer)
 
 ---
 
@@ -149,14 +146,14 @@ A smaller standard deviation indicates that data points are closer to the mean,
 ### `mmean(m0, row)`
 - **Purpose**: Calculate mean along specified columns
 - **Parameters**:
-  - m0: Input matrix (.int)
+  - m0: Input matrix (matrix-handle)
   - column number
 - **Returns**: Mean value (.float)
 
 ### `mstdev(m0, row)`
 - **Purpose**: Calculates standard deviation along specified columns
 - **Parameters**:
-  - m0: Input matrix (.int)
+  - m0: Input matrix (matrix-handle)
   - column number
 - **Returns**: standard deviation value (.float)
 
@@ -168,7 +165,7 @@ stddev = mstdev(id, 2); // Row 2 standard deviation
 ### `colstats`
 - **Purpose**: General statistical analysis
 - **Parameters**:
-  - m0: Input matrix (.int)
+  - m0: Input matrix (matrix-handle)
   - mode: row number 
   - **Returns**: Vector - Matrix(4,1)  with input matrix details
 ```
@@ -339,14 +336,6 @@ mfaplot(loadings, "scree");
    - Memory pool for small matrices
    
 
-3. **Error Handling**
-   ```c
-   #define MATRIX_SUCCESS       0
-   #define MATRIX_MALLOC_ERROR -1
-   #define MATRIX_DIM_ERROR   -2
-   #define MATRIX_SING_ERROR  -3
-   ```
-
 ### Performance Considerations
 1. **Algorithm Selection**
    - Size-based optimization
@@ -373,28 +362,23 @@ mfaplot(loadings, "scree");
   - Convergence verification
   - Error bounds calculation
 
-[Add to Statistical Functions]
-
-
 #### `matrix.stats`
 - **Purpose**: General statistical analysis
 - **Parameters**:
-  - m0: Input matrix (.int)
-  - mode: Analysis mode (.string) - ROW or COL
+  - m0: Input matrix (matrix-handle)
+  - mode: Analysis mode (string) - ROW or COL
 - **Returns**: returns matrix dimension either column or row number 
 
-[Add to Visualization]
-
 #### `matrix.mplot`
 - **Purpose**: Create matrix plot
 - **Parameters**:
-  - m0: Input matrix (.int)
-  - plot_type: Type of plot (.string)
-- **Returns**: Success status (.int)
+  - m0: Input matrix (matrix-handle)
+  - plot_type: Type of plot (string)
+- **Returns**: Success status (integer)
 
 #### `masciiplot`
 - **Purpose**: Create ASCII visualization
 - **Parameters**:
-  - m0: Input matrix (.int)
-  - plot_type: Plot type (.string)
-- **Returns**: Success status (.int)
+  - m0: Input matrix (matrix-handle)
+  - plot_type: Plot type (string)
+- **Returns**: Success status (integer)

lib/plugins/matrix/matrix_test.rexx

@@ -1,113 +1,114 @@
-/* GETPI Plugin Test */
+/* Matrix Plugin Test */
 options levelb
 import rxfnsb
 import matrix
-/*
+
+say '*** Regression Test ***'
 rx=regdat()
 ry=regression(rx,rx+1,'Regression')
-call mprint ry
-say "************* End of Regression *******************"
+call mprint ry,""
+say '*** End of Regression Test ***'
+
+say '*** loading new data, some basic statistical tests ***'
 m1=mcreate(24,5,"Data Matrix")
 call matdata m1
 
-call mprint  m1
+call mprint  m1,""
 say "Mean   COL 1 "mmean(m1,1)
 say "Stddev COL 1 "mstdev(m1,1)
-say "Status M1 "stats(m1,'r')
-say "Status M1 "stats(m1,'c')
+say "Rows of M1 "stats(m1,'r')
+say "Cols of M1 "stats(m1,'c')
 m2=m1
+say '*** end of basic tests ***'
 
+say '*** Create Correlation Matrix ***'
 m9=Mcorr(m1,'Correlation')
-call mprint m9
-say "************* End of Correltation *******************"
-## call mplot m9, "line"
-## call mplot m9, "scatter"
-## call mplot m9, "bar"
-## call mplot m9, "heatmap"
-##  - "line": Line plot
-##  - "scatter": Scatter plot
-##  - "bar": Bar chart
-##  - "heatmap": Heatmap
+call mprint m9,""
+say '*** end of correlation ***'
+
+say '*** Create Covariance Matrix ***'
 ma=Mcov(m1,'Covariance')
-call mprint ma
+call mprint ma,""
+say '*** end of covariance ***'
 
+say '*** Create Correlation Matrix by Matrix operations ***'
 ## lets do the calculation via Matrix operations
 ## 1. transpose Data Matrix
   m2=mtranspose(m1,"Transposed Data Matrix")
-  call mprint m2
+  call mprint m2,""
 say "************* End of Transpose I *******************"
 ## 2. Standardise Data Matrix (mean=0, stddev=1)
   m4=mstandard(m1,"Standardised")
-  call mprint  m4
+  call mprint  m4,""
 ## 3. transpose standardised Data Matrix
   m5=mtranspose(m4,"Transposed Data Matrix")
-  call mprint m5
+  call mprint m5,""
   say "************* End of Transpose II *******************"
 ## 4. multiply Data Matrix with transposed Matrix
   m6=mmult(m5,m4,"close to Correlation Matrix")
-  call mprint m6
+  call mprint m6,""
    say "************* End of Mult II *******************"
 ## 5. Almost there, multiply by rows
   m7=mprod(m6,1/23,"Correlation Matrix")
-  call mprint m7
-    say "************* End of scalar Prod *******************"
-## 6. other stuff: determinant
+  call mprint m7,""
+ say '*** end of correlation ***'
+
+say '*** other functions ***'
    say 'determinante 'mdet(m1)
-    say "************* End of Determinante *******************"
 ## 7. other stuff: L and U Matrix
    mlx=mlu(m6,"L Matrix","U Matrix")
-   call mprint mlx
-   call mprint mlx+1
-stats = mcolstats(m1, "Statistics")
-call mprint(stats)
-
+   call mprint mlx,""
+   call mprint mlx+1,""
+   stats = mcolstats(m1, "Statistics")
+   call mprint stats,""
+say copies('-',72)
+say 'Factor Analysis'
+say copies('-',72)
 say '*** Data Matrix before Factor Analysis ***'
-call mprint m1
-*/
 
 mfdata=mcreate(20,5,'factorial')
-call factdata mfdata
-call mprint mfdata
-mfc=mCorr(mfdata,'Correlation')
-call mprint mfc
+call factdata mfdata    ## load new data
 
+call mprint mfdata,""      ##print it
+mfc=mCorr(mfdata,'Correlation')
+call mprint mfc,'Correlation of loaded data'
+## now do some different facto analysis
 loadings2 = mfactor(mfdata, 2, "Factor Analysis Varimax","Rotate=Varimax,SCORE,diag=3")
-call analyseFact loadings2,mfc
+  call analyseFact loadings2,mfc
 loadings2a = mfactor(mfdata, 1, "Factor Analysis Varimax","Rotate=Varimax,SCORE,diag=3")
-call analyseFact loadings2a,mfc
-
-
+  call analyseFact loadings2a,mfc
 loadings1 = mfactor(mfdata, 2, "Factor Analysis unrotated", "Rotate=none,SCORES,diag=3")
-call analyseFact loadings1,mfc
-
+  call analyseFact loadings1,mfc
 loadings3 = mfactor(mfdata, 2, "Factor Analysis Quartimax","diag=3")
-call analyseFact loadings3,mfc
-
+  call analyseFact loadings3,mfc
 loadings4 = mfactor(mfdata, 2, "Factor Analysis Promax", "Rotate=promax,SCORE,DIAG=5")
-call analyseFact loadings4,mfc
-
+  call analyseFact loadings4,mfc
 call mfree -1            ## free all
 exit
 
 analyseFact: procedure
   arg fact=.int,cor=.int
   say '*** Factors are *** '
-  call mprint fact
+  call mprint fact,""
   factT=mtranspose(fact,"Transposed factorial")
   factmult=Mmult(fact,factT,'Generated Correlation')
-  call mprint factMult
-  call mprint cor
+  call mprint factMult,""
+  call mprint cor,""
+##  call checkCorr cor
   mdif=mcreate(5,5,'Difference FCorrell and Correll')
+
+  rows=stats(cor,'Rows')
+  cols=stats(cor,'cols')
   maxdif=0.0
-  do i=1 to 5
-     do j=1 to 5
+  do i=1 to  rows
+     do j=i+1 to cols
         diff=mget(cor,i,j)-mget(factmult,i,j)
         call mset mdif,i,j,diff
         maxdif=maxdif+abs(diff)
      end
   end
   say "cumulated differences "maxdif
-  call mprint mdif
+##  call mprint mdif,""
 return
 
 
@@ -123,27 +124,54 @@ regression: procedure=.int
 
  ## Step 1: Add a column of ones to X for the intercept
     augmented=mexpand(ind,1,1.0)
-    call mprint augmented
+    call mprint augmented,""
  ## Step 2: calculate ('X*X)
     Xt  = mtranspose(augmented, "Xt");
     XtX = mmult(Xt, augmented,  "XtX");
  ## Step 3: Calculate (X'y)
     XtY = mmult(Xt, dep, "XtY");
-    call mprint xty
+    call mprint xty,""
  ## Step 4: Invert (X'X)
     XtX_inv = minvert(XtX, "XtX_inverted");
     if XtX_inv < 0 then return -21;             ## Singular matrix error
-    call mprint xtx_inv
+    call mprint xtx_inv,""
  ## Step 5: Calculate coefficients: result = (X'X)^(-1)(X'y)
     result = mmult(XtX_inv, XtY, "Regression Coefficients");
-    call mprint result
+    call mprint result,""
  ## Clean up temporary matrices
  say '*************** cleanup *******************'
     say "FREE XT      "xt  mfree(Xt);
     say "FREE XTX     "xtx mfree(XtX);
     say "FREE XTY     "xty mfree(XtY);
     say "FREE XTX_INV "xtx_inv mfree(XtX_inv);
 return result
+
+/* Define thresholds for categorization */
+checkCorr: procedure
+  arg cor=.int
+  strong_threshold = 0.7
+  moderate_threshold = 0.3
+
+  rows=stats(cor,'Rows')
+  cols=stats(cor,'COLs')
+
+/* Interpret the matrix */
+do row = 1 to rows
+    do col = 1 to cols
+        if row = col then iterate
+        corr_value = mget(cor,row,col)
+        if corr_value > strong_threshold then  ,
+           say "Variables "row" and " col "have a strong positive correlation of "corr_value
+        else if corr_value < -strong_threshold then  ,
+           say "Variables "row" and " col "have a strong negative correlation of "corr_value
+        else if corr_value > moderate_threshold then  ,
+           say "Variables "row" and " col "have a moderate positive correlation of "corr_value
+        else if corr_value < -moderate_threshold then  ,
+             say "Variables "row" and " col "have a moderate negative correlation of "corr_value
+    end
+end
+return
+
 factdata:procedure=.int
   arg fact=.int
 call mset fact,1,1,4.2