home *** CD-ROM | disk | FTP | other *** search
/ Fish 'n' More 2 / fishmore-publicdomainlibraryvol.ii1991xetec.iso / fish / applications / xlispstat / src / src2.lzh / XLisp-Stat / matrices2.c < prev    next >
C/C++ Source or Header  |  1990-10-02  |  7KB  |  252 lines
  1. /* matrices2 - Linear algebra routines and Multreg Style Sweep         */
  2. /* Operator. Tolerance should be determined from a keyword variable;   */
  3. /* default tolerance from a global variable.                           */ 
  4. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  5. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  6. /* You may give out copies of this software; for conditions see the    */
  7. /* file COPYING included with this distribution.                       */
  8.  
  9. #include "xlisp.h"
  10. #include "osdef.h"
  11. #ifdef ANSI
  12. #include "xlproto.h"
  13. #include "xlsproto.h"
  14. #else
  15. #include "xlfun.h"
  16. #include "xlsfun.h"
  17. #endif ANSI
  18.  
  19. /* Built in MAKE-SWEEP-MATRIX function */
  20. LVAL xsmakesweepmatrix()
  21. {
  22.   LVAL x, x_data, y, w, x_mean, y_mean, pp, dim, result, result_data;
  23.  
  24.   int n, p, i, j, k;
  25.   double val, dxi, dyi, dv, dw, sum_w, dxik, dxjk, dyj,
  26.   dx_meani, dx_meanj, dy_mean;
  27.   
  28.   /* protect some pointers */
  29.   xlstkcheck(8);
  30.   xlsave(x);
  31.   xlsave(y);
  32.   xlsave(w);
  33.   xlsave(x_mean);
  34.   xlsave(pp);
  35.   xlsave(y_mean);
  36.   xlsave(dim);
  37.   xlsave(result);
  38.   
  39.   x = xsgetmatrix();
  40.   y = xlgetarg();
  41.   w = (moreargs()) ? xlgetarg() : NIL;
  42.   xllastarg();
  43.   
  44.   y = coerce_to_vector(y);
  45.   if (numrows(x) != getsize(y)) xlfail("dimensions do not match");
  46.   if (w != NIL) {
  47.     w = coerce_to_vector(w);
  48.     if (numrows(x) != getsize(w)) xlfail("dimensions do not match");
  49.   }
  50.   x_data = arraydata(x);
  51.  
  52.   /* find dimensions */
  53.   n = getsize(y);
  54.   p = numcols(x);
  55.   
  56.   /* make the result matrix and the mean vector */
  57.   pp = cvfixnum((FIXTYPE) p + 2);
  58.   dim = list2(pp, pp);
  59.   result = newarray(dim, NIL, NIL);
  60.   result_data = arraydata(result);
  61.   x_mean = newvector(p);
  62.   
  63.   /* find the mean of y */
  64.   for (val = 0.0, sum_w = 0.0, i = 0; i < n; i++) {
  65.     dyi = makedouble(getelement(y, i));
  66.     if (w != NIL) {
  67.       dw = makedouble(getelement(w, i));
  68.       sum_w = f_plus(sum_w, dw);
  69.       dyi = f_times(dyi, dw);
  70.     }
  71.     val = f_plus(val, dyi);
  72.   }
  73.   if (w == NIL) sum_w = n;
  74.   if (sum_w <= 0) xlfail("non positive sum of weights");
  75.   y_mean = cvflonum((FLOTYPE) val / sum_w);
  76.   
  77.   /* find the column means */
  78.   for (j = 0; j < p; j++) {
  79.     for (val = 0.0, i = 0; i < n; i++) {
  80.       dxi = makedouble(getelement(x_data, p * i + j));
  81.       if (w != NIL) {
  82.     dw = makedouble(getelement(w, i));
  83.     dxi = f_times(dxi, dw);
  84.       }
  85.       val = f_plus(val, dxi);
  86.     }
  87.     setelement(x_mean, j, cvflonum((FLOTYPE) val / sum_w));
  88.   }
  89.   
  90.   /* put 1/sum_w in topleft, means on left, minus means on top */
  91.   setelement(result_data, 0, cvflonum((FLOTYPE) 1.0 / sum_w));
  92.   for (i = 0; i < p; i++) {
  93.     setelement(result_data, i + 1, 
  94.            cvflonum((FLOTYPE) - getflonum(getelement(x_mean, i))));
  95.     setelement(result_data, (i + 1) * (p + 2), getelement(x_mean, i));
  96.   }
  97.   setelement(result_data, p + 1, cvflonum((FLOTYPE) - getflonum(y_mean)));
  98.   setelement(result_data, (p + 1) * (p + 2), y_mean);
  99.   
  100.   /* put sums of adjusted cross products in body */
  101.   for (i = 0; i < p; i ++) {
  102.     for (j = i; j < p; j++) {
  103.       for (val = 0.0, k = 0; k < n; k++) {
  104.     dxik = makedouble(getelement(x_data, p * k + i));
  105.     dxjk = makedouble(getelement(x_data, p * k + j));
  106.     dx_meani = getflonum(getelement(x_mean, i));
  107.     dx_meanj = getflonum(getelement(x_mean, j));
  108.     dv = f_times(f_minus(dxik, dx_meani), f_minus(dxjk, dx_meanj));
  109.     if (w != NIL) {
  110.       dw = makedouble(getelement(w, k));
  111.       dv = f_times(dv, dw);
  112.     }
  113.     val = f_plus(val, dv);
  114.       }
  115.       setelement(result_data, (i + 1) * (p + 2) + (j + 1), 
  116.          cvflonum((FLOTYPE) val));
  117.       setelement(result_data, (j + 1) * (p + 2) + (i + 1), 
  118.          cvflonum((FLOTYPE) val));
  119.     }
  120.     for (val = 0.0, j = 0; j < n; j++) {
  121.       dxik = makedouble(getelement(x_data, p * j + i));
  122.       dyj = makedouble(getelement(y, j));
  123.       dx_meani = getflonum(getelement(x_mean, i));
  124.       dy_mean = getflonum(y_mean);
  125.       dv = f_times(f_minus(dxik, dx_meani), f_minus(dyj, dy_mean));
  126.       if (w != NIL) {
  127.     dw = makedouble(getelement(w, j));
  128.     dv = f_times(dv, dw);
  129.       }
  130.       val = f_plus(val, dv);
  131.     }
  132.     setelement(result_data, (i + 1) * (p + 2) + (p + 1), 
  133.            cvflonum((FLOTYPE) val));
  134.     setelement(result_data, (p + 1) * (p + 2) + (i + 1), 
  135.            cvflonum((FLOTYPE) val));
  136.   }
  137.   for (val = 0.0, j = 0; j < n; j++) {
  138.     dyj = makedouble(getelement(y, j));
  139.     dy_mean = getflonum(y_mean);
  140.     dv = f_times(f_minus(dyj, dy_mean), f_minus(dyj, dy_mean));
  141.     if (w != NIL) {
  142.       dw = makedouble(getelement(w, j));
  143.       dv = f_times(dv, dw);
  144.     }
  145.     val = f_plus(val, dv);
  146.   }
  147.   setelement(result_data, (p + 1) * (p + 2) + (p + 1), 
  148.          cvflonum((FLOTYPE) val));
  149.   
  150.   /* restore the stack frame */
  151.   xlpopn(8);
  152.   
  153.   return(result);
  154. }
  155.  
  156. /* Internal sweeping algorithm */
  157. int sweep_in_place(a, k, tol)
  158.      LVAL a;
  159.      int k;
  160.      double tol;
  161. {
  162.   LVAL data;
  163.   int rows, cols, i, j;
  164.   double pivot, aij, aik, akj, akk, meps;
  165.   
  166.   if (! matrixp(a)) xlerror("not a matrix", a);
  167.   
  168.   if ((k < 0) || (k >= numrows(a)) || (k >= numcols(a)))
  169.     xlfail("index out of range");
  170.   
  171.   meps = macheps();
  172.   if (tol < meps) tol = meps;
  173.   
  174.   rows = numrows(a);
  175.   cols = numcols(a);
  176.   data = arraydata(a);
  177.  
  178.   pivot = makedouble(getelement(data, cols * k + k));
  179.   
  180.   if ((pivot > tol) || (pivot < -tol)) {
  181.     for (i = 0; i < rows; i++)
  182.       for (j = 0; j < cols; j++)
  183.     if ((i != k) && (j != k)) {
  184.       aij = makedouble(getelement(data, cols * i + j));
  185.       aik = makedouble(getelement(data, cols * i + k));
  186.       akj = makedouble(getelement(data, cols * k + j));
  187.       aij = aij - f_divide(f_times(aik, akj), pivot);
  188.       setelement(data, cols * i + j, cvflonum((FLOTYPE) aij));
  189.     }
  190.     for (i = 0; i < rows; i++) {
  191.       aik = makedouble(getelement(data, cols * i + k));
  192.       if (i != k) {
  193.     aik = f_divide(aik, pivot);
  194.     setelement(data, cols * i + k, cvflonum((FLOTYPE) aik));
  195.       }
  196.     }
  197.     for (j = 0; j < cols; j++) {
  198.       akj = makedouble(getelement(data, cols * k + j));
  199.       if (j != k) {
  200.     akj = -f_divide(akj, pivot);
  201.     setelement(data, cols * k + j, cvflonum((FLOTYPE) akj));
  202.       }
  203.     }
  204.     akk = f_divide(1.0, pivot);
  205.     setelement(data, cols * k + k, cvflonum((FLOTYPE) akk));
  206.     return (1);
  207.   }
  208.   else
  209.     return (0);
  210. }
  211.  
  212. /* Built in SWEEP-OPERATOR function. Should allow in place sweep as option */
  213. LVAL xssweepoperator()
  214. {
  215.   LVAL a, columns, sweptcolumns, tolerances, col;
  216.   LVAL result;
  217.   int c;
  218.   double tol = .000001;
  219.   
  220.   /* create a new stack frame */
  221.   xlstkcheck(3);
  222.   xlsave(sweptcolumns);
  223.   xlsave(tolerances);
  224.   xlsave(result);
  225.   
  226.   a = xsgetmatrix();
  227.   columns = xsgetsequence();
  228.   columns = coerce_to_list(columns);
  229.   tolerances = (moreargs()) ? xsgetsequence() : NIL;
  230.   tolerances = coerce_to_list(tolerances);
  231.   xllastarg();
  232.   
  233.   result = copyarray(a);
  234.   sweptcolumns = NIL;
  235.   
  236.   for (; consp(columns);
  237.        columns = cdr(columns), 
  238.        tolerances = (consp(tolerances)) ? cdr(tolerances) : NIL) {
  239.     col = car(columns);
  240.     if (! fixp(col)) xlerror("not an integer", col);
  241.     c = getfixnum(col);
  242.     tol = (consp(tolerances))  ? makedouble(car(tolerances)) : tol;
  243.     if (sweep_in_place(result, c, tol)) sweptcolumns = cons(col, sweptcolumns);
  244.   }
  245.   result = list2(result, sweptcolumns);
  246.   
  247.   /* restore the previous stack frame */
  248.   xlpopn(3);
  249.   
  250.   return(result);
  251. }
  252.