home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdtrend.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-13  |  25.8 KB  |  781 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdtrend.c    2.15  3/13/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /* grdtrend <input.grd> -N<n_model>[r] [-T<trend.grd>] [-V]
  8.     [-W<weight.grd] [-D<differences.grd]
  9.  
  10. Reads a grdfile and fits a trend surface.  Trend surface
  11. is defined by:
  12.  
  13. m1 +m2*x + m3*y + m4*xy + m5*x*x + m6*y*y + m7*x*x*x
  14.     + m8*x*x*y + m9*x*y*y + m10*y*y*y.
  15.  
  16. n_model is set by the user to be an integer in [1,10]
  17. which sets the number of model coefficients to fit.
  18. Thus:
  19. n_model = 1 gives the mean value of the surface,
  20. n_model = 3 fits a plane,
  21. n_model = 4 fits a bilinear surface,
  22. n_model = 6 fits a biquadratic,
  23. n_model = 10 fits a bicubic surface.
  24.  
  25. The user may write out grdfiles of the fitted surface
  26. [-T<trend.grd>] and / or of the residuals (input data
  27. minus fitted trend) [-D<differences.grd] and / or of
  28. the weights used in iterative fitting [-W<weight.grd].
  29. This last option applies only when the surface is fit
  30. iteratively [-N<n>[r]].
  31.  
  32. A robust fit may be achieved by iterative fitting of
  33. a weighted least squares problem, where the weights
  34. are set according to a scale length based on the 
  35. Median absolute deviation (MAD: Huber, 1982).  The
  36. -N<n>r option acheives this.
  37.  
  38. Author:        W. H. F. Smith
  39. Date:        21 May, 1991.
  40. Version:    1.0, after experience with 1-D case.
  41. Calls:        uses the QR solution of the Normal
  42.         equations furnished by Wm. Menke's
  43.         C routine "gauss".  We gratefully
  44.         acknowledge this contribution.
  45.  
  46. Remarks:
  47.  
  48. We adopt a translation and scaling of the x,y coordinates.
  49. We choose x,y such that they are in [-1,1] over the range
  50. of the grdfile.  If the problem is unweighted, all input
  51. values are filled (no "holes" or NaNs in the input grdfile),
  52. and n_model <= 4 (bilinear or simpler), then the normal
  53. equations matrix (G'G in Menke notation) is diagonal under
  54. this change of coordinates, and the solution is trivial.
  55. In this case, it would be dangerous to try to accumulate
  56. the sums which are the elements of the normal equations;
  57. while they analytically cancel to zero, the addition errors
  58. would likely prevent this.  Therefore we have written a
  59. routine, grd_trivial_model(), to handle this case.
  60.  
  61. If the problem is more complex than the above trivial case,
  62. (missing values, weighted problem, or n_model > 4), then
  63. G'G is not trivial and we just naively accumulate sums in
  64. the G'G matrix.  We hope that the changed coordinates in
  65. [-1,1] will help the accuracy of the problem.  We also use
  66. Legendre polynomials in this case so that the matrix elements
  67. are conveniently sized near 0 or 1.
  68.  
  69. */
  70.  
  71. #include "gmt.h"
  72.  
  73. #define MAX_TABLE_COLS 10    /* Used by Menke routine gauss  */
  74.  
  75. main(argc, argv)
  76. int argc;
  77. char **argv;
  78. {
  79.     int    i, j, k, ierror = 0, iterations, nxy, n_model = 0, dummy[4];
  80.     
  81.     BOOLEAN    error = FALSE, robust = FALSE, trivial, weighted;
  82.  
  83.     double    chisq, old_chisq, zero_test = 1.0e-08, scale = 1.0;
  84.  
  85.     char    *i_filename = NULL, *d_filename = NULL, *t_filename = NULL, *w_filename = NULL;
  86.  
  87.     float    *data;        /* Pointer for array from input grdfile  */
  88.     float    *trend;        /* Pointer for array containing fitted surface  */
  89.     float    *resid;        /* Pointer for array containing residual surface  */
  90.     float    *weight;    /* Pointer for array containing data weights  */
  91.  
  92.     double    *xval;        /* Pointer for array of change of variable:  x[i]  */
  93.     double    *yval;        /* Pointer for array of change of variable:  y[j]  */
  94.     double    *gtg;        /* Pointer for array for matrix G'G normal equations  */
  95.     double    *gtd;        /* Pointer for array for vector G'd normal equations  */
  96.     double    *old;        /* Pointer for array for old model, used for robust sol'n  */
  97.     double    *pstuff;    /* Pointer for array for Legendre polynomials of x[i],y[j]  */
  98.  
  99.     void    set_up_vals();        /* Store x[i], y[j] once for all to save time  */
  100.     void    load_pstuff();        /* Compute Legendre polynomials of x[i],y[j] as needed  */
  101.     void    grd_trivial_model();    /* Fit trivial models.  See Remarks above.  */
  102.     void    compute_trend();    /* Find trend from a model  */
  103.     void    compute_resid();    /* Find residuals from a trend  */
  104.     void    compute_chisq();    /* Find Chi-Squared from weighted residuals  */
  105.     void    compute_robust_weight();    /* Find weights from residuals  */
  106.     void    write_model_parameters();    /* Do reports if gmtdefs.verbose == TRUE  */
  107.     void    load_gtg_and_gtd();        /* Fill normal equations matrices  */
  108.  
  109.     struct GRD_HEADER head_d, head_w;
  110.     
  111. /* Execution begins here with loop over arguments:  */
  112.  
  113.     argc = gmt_begin (argc, argv);
  114.  
  115.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  116.     
  117.     for (i = 1; i < argc; i++) {
  118.         if (argv[i][0] == '-') {
  119.             switch (argv[i][1]) {
  120.  
  121.                 /* Common parameters */
  122.             
  123.                 case 'V':
  124.                 case '\0':
  125.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  126.                     break;
  127.                 
  128.                 /* Supplemental parameters */
  129.             
  130.                 case 'D':
  131.                     d_filename = &argv[i][2];
  132.                     if (d_filename == NULL) {
  133.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -D option:  Must specify file name\n", gmt_program);
  134.                         error = TRUE;
  135.                     }
  136.                     break;
  137.                 case 'N':
  138.                     j = 2;
  139.                     if (argv[i][j] && (argv[i][j] == 'r' || argv[i][j] == 'r') ){
  140.                         robust = TRUE;
  141.                         j++;
  142.                     }
  143.                     if (argv[i][j]) n_model = atoi(&argv[i][j]);
  144.                     break;
  145.                 case 'T':
  146.                     t_filename = &argv[i][2];
  147.                     if (t_filename == NULL) {
  148.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -T option:  Must specify file name\n", gmt_program);
  149.                         error = TRUE;
  150.                     }
  151.                     break;
  152.                 case 'W':
  153.                     w_filename = &argv[i][2];
  154.                     if (w_filename == NULL) {
  155.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -W option:  Must specify file name\n", gmt_program);
  156.                         error = TRUE;
  157.                     }
  158.                     /* OK if this file doesn't exist:  */
  159.                     break;
  160.                 default:
  161.                     error = TRUE;
  162.                     gmt_default_error (argv[i][1]);
  163.                     break;
  164.             }
  165.         }
  166.         else
  167.             i_filename = argv[i];
  168.     }
  169.  
  170.     if (argc == 1 || gmt_quick) {
  171.         fprintf(stderr, "grdtrend %s - Fit trend surface to gridded data\n\n", GMT_VERSION);
  172.         fprintf(stderr,"usage:  grdtrend <input.grd> -N[r]<n_model> [-D<diff.grd>]\n");
  173.         fprintf(stderr,"    [-T<trend.grd>] [-V] [-W<weight.grd>]\n\n");
  174.         
  175.         if (gmt_quick) exit (-1);
  176.         
  177.         fprintf(stderr,"\t<input.grd> = name of grdfile to fit trend to.\n");
  178.         fprintf(stderr,"\t-N[r]<n_model> = # model parameters to fit; integer in [1,10].  Insert r for robust fit.\n");
  179.         fprintf(stderr,"\n\tOPTIONS:\n");
  180.         fprintf(stderr,"\t-D<diff.grd> Supply filename to write grdfile of differences (input - trend).\n");
  181.         fprintf(stderr,"\t-T<trend.grd> Supply filename to write grdfile of trend.\n");
  182.         fprintf(stderr,"\t-V for verbose operation.\n");
  183.         fprintf(stderr,"\t-W<weight.grd> Supply filename if you want to [read and] write grdfile of weights.\n");
  184.         fprintf(stderr,"\t   If <weight.grd> can be read at run, and if robust = FALSE, weighted problem will be solved.\n");
  185.         fprintf(stderr,"\t   If robust = TRUE, weights used for robust fit will be written to <weight.grd>.\n");
  186.         exit(-1);
  187.     }
  188.  
  189.     if (!i_filename) {
  190.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  191.         error++;
  192.     }
  193.     if (n_model <= 0 || n_model > 10) {
  194.         fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Specify 1-10 model parameters\n", gmt_program);
  195.         error++;
  196.     }
  197.  
  198.     if (error) exit (-1);
  199.  
  200. /* End of argument parsing.  */
  201.  
  202.     trivial = ( (n_model < 5) && (!(robust)) && (!w_filename) );
  203.  
  204. /* Read the input file:  */
  205.  
  206.     read_grd_info (i_filename, &head_d);
  207.     grd_init (&head_d, argc, argv, TRUE);
  208.     nxy = head_d.nx * head_d.ny;
  209.     data = (float *) memory (CNULL, nxy, sizeof (float), "grdtrend");
  210.     read_grd (i_filename, &head_d, data, 0.0, 0.0, 0.0, 0.0, dummy, FALSE);
  211.  
  212.     /* Check for NaNs:  */
  213.     i = 0;
  214.     while (trivial && i < nxy) {
  215.         if (bad_float( (double)data[i])) trivial = FALSE;
  216.         i++;
  217.     }
  218.  
  219. /* End input read section.  */
  220.  
  221. /* Allocate other required arrays:  */
  222.  
  223.     trend = (float *) memory (CNULL, nxy, sizeof (float), "grdtrend");
  224.     resid = (float *) memory (CNULL, nxy, sizeof (float), "grdtrend");
  225.     xval = (double *) memory (CNULL, head_d.nx, sizeof (double), "grdtrend");
  226.     yval = (double *) memory (CNULL, head_d.ny, sizeof (double), "grdtrend");
  227.     gtg = (double *) memory (CNULL, n_model*n_model, sizeof (double), "grdtrend");
  228.     gtd = (double *) memory (CNULL, n_model, sizeof (double), "grdtrend");
  229.     old = (double *) memory (CNULL, n_model, sizeof (double), "grdtrend");
  230.     pstuff = (double *) memory (CNULL, n_model, sizeof (double), "grdtrend");
  231.     
  232.  
  233. /* If a weight array is needed, get one:  */
  234.  
  235.     if (weighted = (robust || w_filename) ) {
  236.         weight = (float *) memory (CNULL, nxy, sizeof (float), "grdtrend");
  237.         if (!access (w_filename, R_OK)) {    /* We have weights on input  */
  238.             grd_init (&head_w, argc, argv, FALSE);
  239.             read_grd_info (w_filename, &head_w);
  240.             if (head_w.nx != head_d.nx || head_w.ny != head_d.ny) {
  241.                 fprintf (stderr,"grdtrend:  Input weight file does not match input data file.  Ignoring.\n");
  242.                 for (i = 0; i < nxy; i++) weight[i] = 1.0;
  243.             }
  244.             else
  245.                 read_grd (w_filename, &head_w, weight, 0.0, 0.0, 0.0, 0.0, dummy, FALSE);
  246.         }
  247.         else {
  248.             for (i = 0; i < nxy; i++) weight[i] = 1.0;
  249.         }
  250.     }
  251.  
  252. /* End of weight set up.  */
  253.  
  254. /* Set up xval and yval lookup tables:  */
  255.  
  256.     set_up_vals(xval, head_d.nx, head_d.x_min, head_d.x_max, head_d.x_inc,
  257.          head_d.node_offset);
  258.     set_up_vals(yval, head_d.ny, head_d.y_min, head_d.y_max, head_d.y_inc,
  259.          head_d.node_offset);
  260.  
  261. /* End of set up of lookup values.  */
  262.  
  263. /* Do the problem:  */
  264.  
  265.     if (trivial) {
  266.         grd_trivial_model(data, head_d.nx, head_d.ny, xval, yval, gtd, n_model);
  267.         compute_trend(trend, head_d.nx, head_d.ny, xval, yval, gtd, n_model, pstuff);
  268.         compute_resid(data, trend, resid, nxy);
  269.     }
  270.     else {    /* Problem is not trivial  !!  */
  271.  
  272.         load_gtg_and_gtd(data, head_d.nx, head_d.ny, xval, yval, pstuff, gtg, gtd, n_model, weight, weighted);
  273.         gauss(gtg, gtd, n_model, n_model, zero_test, &ierror, 1);
  274.         if (ierror) {
  275.             fprintf(stderr,"grdtrend:  Gauss returns error code %d\n", ierror);
  276.             exit(-1);
  277.         }
  278.         compute_trend(trend, head_d.nx, head_d.ny, xval, yval, gtd, n_model, pstuff);
  279.         compute_resid(data, trend, resid, nxy);
  280.  
  281.         if (robust) {
  282.             compute_chisq(resid, weight, nxy, &chisq, scale);
  283.             iterations = 1;
  284.             do {
  285.                 old_chisq = chisq;
  286.                 for(k = 0; k < n_model; k++) old[k] = gtd[k];
  287.                 compute_robust_weight(resid, weight, nxy, &scale);
  288.                 load_gtg_and_gtd(data, head_d.nx, head_d.ny, xval, yval, pstuff, gtg, gtd, n_model, weight, weighted);
  289.                 gauss(gtg, gtd, n_model, n_model, zero_test, &ierror, 1);
  290.                 if (ierror) {
  291.                     fprintf(stderr,"grdtrend:  Gauss returns error code %d\n", ierror);
  292.                     exit(-1);
  293.                 }
  294.                 compute_trend(trend, head_d.nx, head_d.ny, xval, yval, gtd, n_model, pstuff);
  295.                 compute_resid(data, trend, resid, nxy);
  296.                 compute_chisq(resid, weight, nxy, &chisq, scale);
  297.                 if (gmtdefs.verbose) fprintf(stderr,"grdtrend Robust iteration %d:  Old Chi Squared:  %.8lg  New Chi Squared %.8lg\n", iterations, old_chisq, chisq);
  298.                 iterations++;
  299.             } while ( (old_chisq / chisq) > 1.0001);
  300.  
  301.             /* Get here when new model not significantly better; use old one:  */
  302.  
  303.             for(k = 0; k < n_model; k++) gtd[k] = old[k];
  304.             compute_trend(trend, head_d.nx, head_d.ny, xval, yval, gtd, n_model, pstuff);
  305.             compute_resid(data, trend, resid, nxy);
  306.         }
  307.     }
  308.  
  309. /* End of do the problem section.  */
  310.  
  311. /* Get here when ready to do output:  */
  312.     
  313.     if (gmtdefs.verbose) write_model_parameters(gtd, n_model);
  314.     if (t_filename) {
  315.         strcpy (head_d.title, "trend surface");
  316.         write_grd (t_filename, &head_d, trend, 0.0, 0.0, 0.0, 0.0, dummy, FALSE);
  317.     }
  318.     if (d_filename) {
  319.         strcpy (head_d.title, "trend residuals");
  320.         write_grd (d_filename, &head_d, resid, 0.0, 0.0, 0.0, 0.0, dummy, FALSE);
  321.     }
  322.     if (w_filename) {
  323.         strcpy (head_d.title, "trend weights");
  324.         write_grd (d_filename, &head_d, weight, 0.0, 0.0, 0.0, 0.0, dummy, FALSE);
  325.     }
  326.  
  327. /* That's all, folks!  */
  328.  
  329.     if (weighted) free((char *)weight);
  330.     free((char *)pstuff);
  331.     free((char *)gtd);
  332.     free((char *)gtg);
  333.     free((char *)yval);
  334.     free((char *)xval);
  335.     free((char *)resid);
  336.     free((char *)trend);
  337.     free((char *)data);
  338.  
  339.     gmt_end (argc, argv);
  340. }        
  341.  
  342. void  set_up_vals(val, nval, vmin, vmax, dv, pixel_reg)
  343. double    val[], vmin, vmax, dv;
  344. int    nval, pixel_reg;
  345. {
  346.     int    i;
  347.     double  v, middle, drange, true_min, true_max;
  348.  
  349.     true_min = (pixel_reg) ? vmin + 0.5 * dv : vmin;
  350.     true_max = (pixel_reg) ? vmax - 0.5 * dv : vmax;
  351.  
  352.     middle = 0.5 * (true_min + true_max);
  353.     drange = 2.0 / (true_max - true_min);
  354.     for (i = 0; i < nval; i++) {
  355.         v = true_min + i * dv;
  356.         val[i] = (v - middle) * drange;
  357.     }
  358.     /* Just to be sure no rounding outside:  */
  359.     val[0] = -1.0;
  360.     val[nval - 1] = 1.0;
  361.     return;
  362. }
  363.  
  364. void    load_pstuff(pstuff, n_model, x, y, newx, newy)
  365. double    pstuff[], x, y;
  366. int    n_model, newx, newy;
  367. {
  368.     /* If either x or y has changed, compute new Legendre polynomials as needed  */
  369.  
  370.     if (newx) {
  371.         if (n_model >= 2) pstuff[1] = x;
  372.         if (n_model >= 5) pstuff[4] = 0.5*(3.0*pstuff[1]*pstuff[1] - 1.0);
  373.         if (n_model >= 7) pstuff[6] = (5.0*pstuff[1]*pstuff[4] - 2.0*pstuff[1])/3.0;
  374.     }
  375.     if (newy) {
  376.         if (n_model >= 3) pstuff[2] = y;
  377.         if (n_model >= 6) pstuff[5] = 0.5*(3.0*pstuff[2]*pstuff[2] - 1.0);
  378.         if (n_model >= 10) pstuff[9] = (5.0*pstuff[2]*pstuff[5] - 2.0*pstuff[2])/3.0;
  379.     }
  380.     /* In either case, refresh cross terms:  */
  381.  
  382.     if (n_model >= 4) pstuff[3] = pstuff[1]*pstuff[2];
  383.     if (n_model >= 8) pstuff[7] = pstuff[4]*pstuff[2];
  384.     if (n_model >= 9) pstuff[8] = pstuff[1]*pstuff[5];
  385.  
  386.     return;
  387. }
  388.  
  389. void    compute_trend(trend, nx, ny, xval, yval, gtd, n_model, pstuff)
  390. float    trend[];
  391. double    xval[], yval[], gtd[], pstuff[];
  392. int    nx, ny, n_model;
  393. {
  394.     int    i, j, k, ij;
  395.  
  396.     for (ij = 0, j = 0; j < ny; j++) {
  397.         for (i = 0; i < nx; i++, ij++) {
  398.             load_pstuff(pstuff, n_model, xval[i], yval[j], 1, (!(i)));
  399.             trend[ij] = gtd[0];
  400.             for (k = 1; k < n_model; k++) {
  401.                 trend[ij] += (pstuff[k]*gtd[k]);
  402.             }
  403.         }
  404.     }
  405. }
  406.  
  407. void    compute_resid(data, trend, resid, nxy)
  408. float    data[], trend[], resid[];
  409. int    nxy;
  410. {
  411.     int    i;
  412.  
  413.     for (i = 0; i < nxy; i++) {
  414.         if (bad_float( (double)data[i])) {
  415.             resid[i] = data[i];
  416.         }
  417.         else {
  418.             resid[i] = data[i] - trend[i];
  419.         }
  420.     }
  421.     return;
  422. }
  423.  
  424. void    grd_trivial_model(data, nx, ny, xval, yval, gtd, n_model)
  425. float    data[];
  426. double    xval[], yval[], gtd[];
  427. int    nx, ny, n_model;
  428. {
  429.     /* Routine to fit up elementary polynomial model of grd data, 
  430.     model = gtd[0] + gtd[1]*x + gtd[2]*y + gtd[3] * x * y,
  431.     where x,y are normalized to range [-1,1] and there are no
  432.     NaNs in grdfile, and problem is unweighted least squares.  */
  433.  
  434.     int    i, j, ij;
  435.     double    x2, y2, sumx2 = 0.0, sumy2 = 0.0, sumx2y2 = 0.0;
  436.  
  437.     /* First zero the model parameters to use for sums:  */
  438.  
  439.     for (i = 0; i < n_model; i++) gtd[i] = 0.0;
  440.  
  441.     /* Now accumulate sums:  */
  442.  
  443.     for (ij = 0, j = 0; j < ny; j++) {
  444.         y2 = yval[j] * yval[j];
  445.         for (i = 0; i < nx; i++, ij++) {
  446.             x2 = xval[i] * xval[i];
  447.             sumx2 += x2;
  448.             sumy2 += y2;
  449.             sumx2y2 += (x2 * y2);
  450.             gtd[0] += data[ij];
  451.             if (n_model >= 2) gtd[1] += data[ij] * xval[i];
  452.             if (n_model >= 3) gtd[2] += data[ij] * yval[j];
  453.             if (n_model == 4) gtd[3] += data[ij] * xval[i] * yval[j];
  454.         }
  455.     }
  456.  
  457.     /* See how trivial it is?  */
  458.  
  459.     gtd[0] /= (nx * ny);
  460.     if (n_model >= 2) gtd[1] /= sumx2;
  461.     if (n_model >= 3) gtd[2] /= sumy2;
  462.     if (n_model == 4) gtd[3] /= sumx2y2;
  463.  
  464.     return;
  465. }
  466.  
  467. void    compute_chisq(resid, weight, nxy, chisq, scale)
  468. float    resid[], weight[];
  469. int    nxy;
  470. double    *chisq, scale;
  471. {
  472.     int    i;
  473.     double    tmp;
  474.  
  475.     *chisq = 0.0;
  476.     for (i = 0; i < nxy; i++) {
  477.         if (bad_float( (double)resid[i]))continue;
  478.         tmp = resid[i];
  479.         if (scale != 1.0) tmp /= scale;
  480.         tmp *= tmp;
  481.  
  482.         if (weight[i] == 1.0) {
  483.             *chisq += tmp;
  484.         }
  485.         else {
  486.             /* Weight has already been squared  */
  487.             *chisq += (tmp * weight[i]);
  488.         }
  489.     }
  490.     return;
  491. }
  492.  
  493. void    compute_robust_weight(resid, weight, nxy, scale)
  494. float    resid[], weight[];
  495. int    nxy;
  496. double    *scale;
  497. {
  498.     int    i, j;
  499.     double    r, mad;
  500.  
  501.     for (i = j = 0; i < nxy; i++) {
  502.         if (bad_float( (double)resid[i]))continue;
  503.         weight[j] = fabs((double)resid[i]);
  504.         j++;
  505.     }
  506.  
  507.     qsort( (char *)weight, j, sizeof(float), comp_float_asc);
  508.  
  509.     if (j%2) {
  510.         mad = weight[j/2];
  511.     }
  512.     else {
  513.         mad = 0.5 *(weight[j/2] + weight[j/2 - 1]);
  514.     }
  515.  
  516.     /* Adjust mad to equal Gaussian sigma:  */
  517.  
  518.     *scale = 1.4826 * mad;
  519.  
  520.     /* Use weight according to Huber (1981), but squared:  */
  521.  
  522.     for (i = 0; i < nxy; i++) {
  523.         if (bad_float( (double)resid[i])) {
  524.             weight[i] = resid[i];
  525.             continue;
  526.         }
  527.         r = fabs(resid[i]) / (*scale);
  528.         
  529.         weight[i] = (r <= 1.5) ? 1.0 : (3.0 - 2.25/r) / r;
  530.     }
  531.     return;
  532. }
  533.  
  534. void    write_model_parameters(gtd, n_model)
  535. double    gtd[];
  536. int    n_model;
  537. {
  538.     int    i;
  539.     char    pbasis[10][12];
  540.  
  541.     sprintf(pbasis[0], "Mean\0");
  542.     sprintf(pbasis[1], "X\0");
  543.     sprintf(pbasis[2], "Y\0");
  544.     sprintf(pbasis[3], "X*Y\0");
  545.     sprintf(pbasis[4], "P2(x)\0");
  546.     sprintf(pbasis[5], "P2(y)\0");
  547.     sprintf(pbasis[6], "P3(x)\0");
  548.     sprintf(pbasis[7], "P2(x)*P1(y)\0");
  549.     sprintf(pbasis[8], "P1(x)*P2(y)\0");
  550.     sprintf(pbasis[9], "P3(y)\0");
  551.  
  552.     for(i = 0; i < n_model; i++) {
  553.         fprintf(stderr,"Coefficient fit to %s:  %.8lg\n", pbasis[i], gtd[i]);
  554.     }
  555.     return;
  556. }
  557.  
  558. void    load_gtg_and_gtd(data, nx, ny, xval, yval, pstuff, gtg, gtd, n_model, weight, weighted)
  559. float    data[], weight[];
  560. int    nx, ny, n_model, weighted;
  561. double    xval[], yval[], pstuff[], gtg[], gtd[];
  562. {
  563.     /* Routine to load the matrix G'G (gtg) and vector G'd (gtd)
  564.     for the normal equations.  Routine uses indices i,j to refer
  565.     to the grdfile of data, and k,l to refer to the k_row, l_col
  566.     of the normal equations matrix.  We need sums of [weighted]
  567.     data and model functions in gtg and gtd.  We save time by
  568.     loading only lower triangular part of gtg and then filling
  569.     by symmetry after i,j loop.  */
  570.  
  571.     int    i, j, ij, k, l, n_used;
  572.  
  573. /*    First zero things out to start:  */
  574.  
  575.     n_used = 0;
  576.     for (k = 0; k < n_model; k++) {
  577.         gtd[k] = 0.0;
  578.         for (l = 0; l < n_model; l++) {
  579.             gtg[k*n_model+l] = 0.0;
  580.         }
  581.     }
  582.  
  583. /*  Now get going.  Have to load_pstuff separately in i and j,
  584.     because it is possible that we skip data when i = 0.
  585.     Loop over all data:  */
  586.  
  587.     for (ij = 0, j = 0; j < ny; j++ ) {
  588.         load_pstuff(pstuff, n_model, xval[0], yval[j], 0, 1);
  589.         for (i = 0; i < nx; i++, ij++) {
  590.  
  591.             if (bad_float( (double)data[ij]))continue;
  592.  
  593.             n_used++;
  594.             load_pstuff(pstuff, n_model, xval[i], yval[j], 1, 0);
  595.  
  596. /* If weighted  */    if (weighted) {
  597.                 /* Loop over all gtg and gtd elements:  */
  598.                 gtd[0] += (data[ij] * weight[ij]);
  599.                 gtg[0] += (weight[ij]);
  600.                 for (k = 1; k < n_model; k++) {
  601.                     gtd[k] += (data[ij] * weight[ij] * pstuff[k]);
  602.                     gtg[k] += (weight[ij] * pstuff[k]);
  603.                     for (l = k; l < n_model; l++) {
  604.                         gtg[k + l*n_model] += (pstuff[k]*pstuff[l]*weight[ij]);
  605.                     }
  606.                 }
  607.             }
  608. /* If !weighted  */    else {
  609.                 /* Loop over all gtg and gtd elements:  */
  610.                 gtd[0] += data[ij];
  611.                 for (k = 1; k < n_model; k++) {
  612.                     gtd[k] += (data[ij] * pstuff[k]);
  613.                     gtg[k] += pstuff[k];
  614.                     for (l = k; l < n_model; l++) {
  615.                         gtg[k + l*n_model] += (pstuff[k]*pstuff[l]);
  616.                     }
  617.                 }
  618. /* End if  */        }
  619.         }
  620.     }
  621. /* End of loop over data i,j  */
  622.  
  623. /* Now if !weighted, use more accurate sum for gtg[0], and set symmetry:  */
  624.  
  625.     if (!(weighted)) gtg[0] = n_used;
  626.  
  627.     for (k = 0; k < n_model; k++) {
  628.         for (l = 0; l < k; l++) {
  629.             gtg[l + k*n_model] = gtg[k + l*n_model];
  630.         }
  631.     }
  632. /* That is all there is to it!  */
  633.  
  634.     return;
  635. }
  636.                 
  637. gauss(a,vec,n,nstore,test,ierror,itriag)
  638. double *a, vec[], test;
  639. int n, nstore, *ierror, itriag;
  640. {
  641.  
  642. /* subroutine gauss, by william menke */
  643. /* july 1978 (modified feb 1983, nov 85) */
  644.  
  645. /* a subroutine to solve a system of n linear equations in n unknowns*/
  646. /* where n doesn't exceed MAX_TABLE_COLS */
  647. /* gaussian reduction with partial pivoting is used */
  648. /*      a               (sent, destroyed)       n by n matrix           */
  649. /*      vec             (sent, overwritten)     n vector, replaced w/ solution*/
  650. /*      nstore          (sent)                  dimension of a  */
  651. /*      test            (sent)                  div by zero check number*/
  652. /*      ierror          (returned)              zero on no error*/
  653. /*      itriag          (sent)                  matrix triangularized only*/
  654. /*                                               on TRUE useful when solving*/
  655. /*                                               multiple systems with same a */
  656.         static int isub[MAX_TABLE_COLS], l1;
  657.         int line[MAX_TABLE_COLS], iet, ieb, i, j, k, l, j2;
  658.         double big, testa, b, sum;
  659.         
  660.  
  661.         iet=0;  /* initial error flags, one for triagularization*/
  662.         ieb=0;  /* one for backsolving */
  663.  
  664. /* triangularize the matrix a*/
  665. /* replacing the zero elements of the triangularized matrix */
  666. /* with the coefficients needed to transform the vector vec */
  667.  
  668.         if (itriag) {   /* triangularize matrix */
  669.  
  670.                 for( j=0; j<n; j++ ) {      /*line is an array of flags*/
  671.                         line[j]=0; 
  672.                         /* elements of a are not moved during pivoting*/
  673.                         /* line=0 flags unused lines */
  674.                         }    /*end for j*/
  675.                         
  676.                 for( j=0; j<n-1; j++ ) {
  677.                         /*  triangularize matrix by partial pivoting */
  678.                        big = 0.0; /* find biggest element in j-th column*/
  679.                                   /* of unused portion of matrix*/
  680.                        for( l1=0; l1<n; l1++ ) {
  681.                                if( line[l1]==0 ) {
  682.                                        testa=(double) fabs(
  683.                                                 (double) (*(a+l1*nstore+j)) );
  684.                                        if (testa>big) {
  685.                                                 i=l1;
  686.                                                 big=testa;
  687.                                                 } /*end if*/
  688.                                         } /*end if*/
  689.                                 } /*end for l1*/
  690.                        if( big<=test) {   /* test for div by 0 */
  691.                                iet=1;
  692.                                } /*end if*/
  693.  
  694.                        line[i]=1;  /* selected unused line becomes used line */
  695.                        isub[j]=i;  /* isub points to j-th row of tri. matrix */
  696.  
  697.                        sum=1.0/(*(a+i*nstore+j)); 
  698.                                 /*reduce matrix towards triangle */
  699.                        for( k=0; k<n; k++ ) {
  700.                                 if( line[k]==0 ) {
  701.                                         b=(*(a+k*nstore+j))*sum;
  702.                                         for( l=j+1; l<n; l++ ) {
  703.                                                *(a+k*nstore+l)=
  704.                                                         (*(a+k*nstore+l))
  705.                                                         -b*(*(a+i*nstore+l));
  706.                                                } /*end for l*/
  707.                                        *(a+k*nstore+j)=b;
  708.                                         } /*end if*/
  709.                                 } /*end for k*/
  710.                         } /*end for j*/
  711.  
  712.                for( j=0; j<n; j++ ) {
  713.                         /*find last unused row and set its pointer*/
  714.                         /*  this row contians the apex of the triangle*/
  715.                         if( line[j]==0) {
  716.                                 l1=j;   /*apex of triangle*/
  717.                                 isub[n-1]=j;
  718.                                 break;
  719.                                 } /*end if*/
  720.                         } /*end for j*/
  721.  
  722.                 } /*end if itriag true*/
  723.                 
  724.         /*start backsolving*/
  725.         
  726.         for( i=0; i<n; i++ ) {  /* invert pointers. line(i) now gives*/
  727.                                 /* row no in triang matrix of i-th row*/
  728.                                 /* of actual matrix */
  729.                 line[isub[i]] = i;
  730.                 } /*end for i*/
  731.  
  732.         for( j=0; j<n-1; j++) { /*transform the vector to match triang. matrix*/
  733.                b=vec[isub[j]];
  734.                for( k=0; k<n; k++ ) {
  735.                       if (line[k]>j) {  /* skip elements outside of triangle*/
  736.                                 vec[k]=vec[k]-(*(a+k*nstore+j))*b;
  737.                                 } /*end if*/
  738.                         } /*end for k*/
  739.                 } /*end for j*/
  740.  
  741.       b = *(a+l1*nstore+(n-1));   /*apex of triangle*/
  742.       if( ((double)fabs( (double) b))<=test) {
  743.                 /*check for div by zero in backsolving*/
  744.                 ieb=2;
  745.                 } /*end if*/
  746.       vec[isub[n-1]]=vec[isub[n-1]]/b;
  747.  
  748.       for( j=n-2; j>=0; j-- ) { /* backsolve rest of triangle*/
  749.                 sum=vec[isub[j]];
  750.                 for( j2=j+1; j2<n; j2++ ) {
  751.                         sum = sum - vec[isub[j2]] * (*(a+isub[j]*nstore+j2));
  752.                         } /*end for j2*/
  753.                         b = *(a+isub[j]*nstore+j);
  754.                if( ((double)fabs((double)b))<=test) {
  755.                         /* test for div by 0 in backsolving */
  756.                         ieb=2;
  757.                         } /*end if*/
  758.                 vec[isub[j]]=sum/b;   /*solution returned in vec*/
  759.                 } /*end for j*/
  760.  
  761. /*put the solution vector into the proper order*/
  762.  
  763.       for( i=0; i<n; i++ ) {    /* reorder solution */
  764.                 for( k=i; k<n; k++ ) {  /* search for i-th solution element */
  765.                         if( line[k]==i ) {
  766.                                 j=k;
  767.                                 break;
  768.                                 } /*end if*/
  769.                         } /*end for k*/
  770.                b = vec[j];       /* swap solution and pointer elements*/
  771.                vec[j] = vec[i];
  772.                vec[i] = b;
  773.                line[j] = line[i];
  774.                 } /*end for i*/
  775.  
  776.       *ierror = iet + ieb;   /* set final error flag*/
  777. }
  778.  
  779.  
  780.  
  781.