home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / trend2d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-13  |  24.7 KB  |  740 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)trend2d.c    2.12  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. /*
  8.  * trend2d.c [<xyz[w]file>] -F<output_flags> -N<n_m_parameters>[r] 
  9.  *    [-C<condition_#>] [-I[<confid>]] [-V] [-W]
  10.  *
  11.  * where:
  12.  *    [<xyz[w]file>] is an ascii file with x y z in first 3 columns
  13.  *        [or x y z w in first 4 columns].  Default reads from stdin.
  14.  *    -F<output_flags> is a string of at least one, up to five, in
  15.  *        and order, from the set {x y z m r w}.  x,y,z = input,
  16.  *        m = model, r = residual = z-m, and w= weight used.
  17.  *    -N<n_m_parameters>[r]
  18.  *        If iterative Robust fitting desired, use -N<#>r, else -N.
  19.  *        [Max] Number of terms in the model is <n_m_parameters>.
  20.  *        Example:  Robust bilinear surface:  -N4r.  Max n = 10.
  21.  *    [-C<condition_#>] Cut off eigenvalue spectrum; use only eigen-
  22.  *        values such that (lambda_max / lambda[i]) < condition_#.
  23.  *    [-I[<confid>]] Iteratively Increment the number of model parameters,
  24.  *        searching for the significant model size, up to a maximum
  25.  *        set by <n_m_parameters>.  We start with a 1 parameter
  26.  *        model and then iteratively increase the number of
  27.  *        model parameters, m, while m <= <n_m_parameters> &&
  28.  *        reduction in variance from i to i+1 is significant
  29.  *        at the <confid> level according to F test.  If user sets
  30.  *        -I without giving <confid> then <confid> = 0.95.
  31.  *    [-V]    Verbose operation.
  32.  *    [-W]    Weighted data are input.  Read 4 cols and use 4th as weight.
  33.  *
  34.  *
  35.  * Read stdin or file of x y z triples, or weighted data, x y z w.  Fit 
  36.  * a regression model z = f(x,y) + e, where e are error misfits and f(x,y)
  37.  * has some user-prescribed functional form.  The user may choose the number
  38.  * of terms in the model to fit, whether to seek iterative refinement robust
  39.  * w.r.t. outliers, and whether to seek automatic discovery of the significant
  40.  * number of model parameters.
  41.  *
  42.  * Adapted from trend1d by w. h. f. smith.
  43.  *
  44.  *
  45.  * During model fitting the data x,y coordinates are normalized into the domain
  46.  * [-1, 1] for Chebyshev Polynomial fitting.  Before writing out the data the
  47.  * coordinates are rescaled to match the original input values.
  48.  *
  49.  *
  50.  * Author:    W. H. F. Smith
  51.  * Date:    17 June 1991-1995.
  52.  */
  53.  
  54. #include "gmt.h"
  55.  
  56. #define N_OUTPUT_CHOICES 6
  57.  
  58. struct    DATA {
  59.     double    x;
  60.     double    y;
  61.     double    z;
  62.     double    m;
  63.     double    r;
  64.     double    w;
  65. }    *data;
  66.  
  67. char format[20];
  68.  
  69. main(argc,argv)
  70. int    argc;
  71. char    **argv;
  72. {
  73.     void    read_data(), write_output(), transform_x(), untransform_x(), recompute_weights();
  74.     void    allocate_array_space(), free_the_memory(), calc_m_and_r(), move_model_a_to_b();
  75.     void    load_gtg_and_gtd(), solve_system();
  76.  
  77.     int    i, j, n_data, n_outputs, n_model, n_model_max, np, significant, rank;
  78.     BOOLEAN    error = FALSE, weighted_input = FALSE, weighted_output = FALSE, robust = FALSE, increment = FALSE;
  79.  
  80.     double    c_no = 1.0e06;    /* Condition number for matrix solution  */
  81.     double    confid = 0.51;    /* Confidence interval for significance test  */
  82.     double    *gtg, *v, *gtd, *lambda, *workb, *workz, *c_model, *o_model, *w_model, *work;    /* Arrays  */
  83.     double    xmin, xmax, ymin, ymax, c_chisq, o_chisq, w_chisq, scale = 1.0, prob;
  84.     double    get_chisq();
  85.  
  86.     char    output_choice[N_OUTPUT_CHOICES];
  87.     FILE    *fp = NULL;
  88.  
  89.     argc = gmt_begin (argc, argv);
  90.  
  91.     n_outputs = 0;
  92.     n_model_max = 0;
  93.         for (i = 0; i < N_OUTPUT_CHOICES; i++) output_choice[i] = 0;
  94.     sprintf (format, "%s\t\0", gmtdefs.d_format);
  95.  
  96.         for (i = 1; i < argc; i++) {
  97.                 if (argv[i][0] == '-') {
  98.                         switch (argv[i][1]) {
  99.                 /* Common parameters */
  100.             
  101.                 case 'H':
  102.                 case 'V':
  103.                 case ':':
  104.                 case '\0':
  105.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  106.                     break;
  107.                 
  108.                 /* Supplemental parameters */
  109.                 
  110.  
  111.                 case 'F':
  112.                     j = 2;
  113.                                        while (argv[i][j]) {
  114.                                                 switch (argv[i][j]) {
  115.                                                         case 'x':
  116.                                                                 output_choice[j-2] = 'x';
  117.                                                                 break;
  118.                                                         case 'y':
  119.                                                                 output_choice[j-2] = 'y';
  120.                                                                 break;
  121.                                                         case 'z':
  122.                                                                 output_choice[j-2] = 'z';
  123.                                                                 break;
  124.                                                         case 'm':
  125.                                                                 output_choice[j-2] = 'm';
  126.                                                                 break;
  127.                                                         case 'r':
  128.                                                                 output_choice[j-2] = 'r';
  129.                                                                 break;
  130.                                                         case 'w':
  131.                                                                 output_choice[j-2] = 'w';
  132.                                 weighted_output = TRUE;
  133.                                                                 break;
  134.                                                         default:
  135.                                                                 error = TRUE;
  136.                                  fprintf (stderr, "%s: GMT SYNTAX ERROR -F option.  Unrecognized output choice %c\n", gmt_program, argv[i][j]);
  137.                                  break;
  138.                                                 }
  139.                                                 n_outputs++;
  140.                                                 j++;
  141.                                         }
  142.                     break;
  143.                 case 'C':
  144.                     c_no = atof(&argv[i][2]);
  145.                     break;
  146.                 case 'I':
  147.                     increment = TRUE;
  148.                     confid = (argv[i][2]) ? atof(&argv[i][2]) : 0.51;
  149.                     break;
  150.                 case 'N':
  151.                     if (argv[i][strlen (argv[i]) - 1] == 'r') robust = TRUE;
  152.                     n_model_max = (argv[i][2]) ? atoi(&argv[i][2]) : 0;
  153.                     break;
  154.                 case 'W':
  155.                     weighted_input = TRUE;
  156.                     break;
  157.                                 default:
  158.                                         error = TRUE;
  159.                     gmt_default_error (argv[i][1]);
  160.                                         break;
  161.                         }
  162.                 }
  163.                 else {
  164.                         if ((fp = fopen(argv[i], "r")) == NULL) {
  165.                                 fprintf (stderr, "trend2d:  Could not open file %s\n", argv[i]);
  166.                                 error = TRUE;
  167.                         }
  168.                 }
  169.         }
  170.  
  171.     if (argc == 1 || gmt_quick) {
  172.         fprintf(stderr,"trend2d %s - Fit a [weighted] [robust] polynomial for z = f(x,y) to ascii xyz[w]\n\n", GMT_VERSION);
  173.         fprintf(stderr,"usage:  trend2d -N<n_model>[r] [<xyz[w]file>] [-F<xymrw>] [-C<condition_#>] [-I[<confidence>]] [-V] [-W]\n\n");
  174.         
  175.         if (gmt_quick) exit (-1);
  176.  
  177.         fprintf(stderr,"\t-N<n_model>[r] fit a [robust] model with <n_model> terms.  <n_model> in [1,10].  E.g., robust quadratic = -N3r.\n");
  178.         fprintf (stderr, "\n\tOPTIONS:\n");
  179.         fprintf(stderr,"\t[<xyz[w]file>] name of ascii file, first 2 cols = x y [3 cols = x y w].  [Default reads stdin].\n");
  180.         fprintf(stderr,"\t-F<xymrw>  Choose at least 1, up to 6, any order, of xyzmrw for ascii output to stdout.\n");
  181.         fprintf(stderr,"\t\tx=x, y=y, z=z, m=model, r=residual=z-m, w=weight.  w determined iteratively if robust fit used.\n");
  182.         fprintf(stderr,"\t-C  Truncate eigenvalue spectrum so matrix has <condition_#>.  [Default = 1.0e06].\n");
  183.         fprintf(stderr,"\t-I  Iteratively Increase # model parameters, to a max of <n_model> so long as the\n");
  184.         fprintf(stderr,"\t\treduction in variance is significant at the <confidence> level.\n");
  185.         fprintf(stderr,"\t\tGive -I without a number to default to 0.51 confidence level.\n");
  186.         fprintf(stderr,"\t-V  Verbose operation.  Program reports progress to stderr.\n");
  187.         fprintf(stderr,"\t-W  Weighted input given, weights in 4th column.  [Default is unweighted].\n\n");
  188.         exit(-1);
  189.     }
  190.  
  191.     if (c_no <= 1.0) {
  192.         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option.  Condition number must be larger than unity\n", gmt_program);
  193.         error++;
  194.     }    
  195.     if (confid < 0.0 || confid > 1.0) {
  196.         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option.  Give 0 < confidence level < 1.0\n", gmt_program);
  197.         error++;
  198.     }    
  199.     if (n_outputs > N_OUTPUT_CHOICES) {
  200.         fprintf (stderr, "%s: GMT SYNTAX ERROR -F option.  Too many output columns specified (%d)\n", gmt_program, n_outputs);
  201.         error++;
  202.     }    
  203.     if (n_model_max <= 0 || n_model_max > 10) {
  204.         fprintf (stderr, "%s: GMT SYNTAX ERROR -N option.  Must request 1-10 parameters\n", gmt_program);
  205.         error++;
  206.     }    
  207.  
  208.     if (error) exit (-1);
  209.  
  210.     np = n_model_max;    /* Row dimension for matrices gtg and v  */
  211.     allocate_array_space(np, >g, &v, >d, &lambda, &workb, &workz, &c_model, &o_model, &w_model);
  212.  
  213.     read_data(&data, &n_data, &xmin, &xmax, &ymin, &ymax, weighted_input, &work, fp);
  214.  
  215.     if (xmin == xmax || ymin == ymax) {
  216.         fprintf(stderr,"trend2d:  Fatal error in input data.  X min = X max.\n");
  217.         exit(-1);
  218.     }
  219.     if (n_data == 0) {
  220.         fprintf(stderr,"trend2d:  Fatal error.  Could not read any data.\n");
  221.         exit(-1);
  222.     }
  223.     if (n_data < n_model_max) {
  224.         fprintf(stderr,"trend2d:  Warning.  Ill-posed problem.  n_data < n_model_max.\n");
  225.     }
  226.  
  227.     transform_x(data, n_data, xmin, xmax, ymin, ymax);    /* Set domain to [-1, 1] or [-pi, pi]  */
  228.  
  229.     if (gmtdefs.verbose) {
  230.         fprintf(stderr,"trend2d:  Read %d data with X values from %.8lg to %.8lg\n", n_data, xmin, xmax);
  231.         fprintf(stderr,"N_model\tRank\tChi_Squared\tSignificance\n");
  232.     }
  233.  
  234.     if (increment) {
  235.         n_model = 1;
  236.  
  237.         /* Fit first model  */
  238.         load_gtg_and_gtd(data, n_data, gtg, gtd, workb, n_model, np);
  239.         solve_system(gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, c_no, &rank);
  240.         calc_m_and_r(data, n_data, c_model, n_model, workb);
  241.         c_chisq = get_chisq(data, n_data, n_model);
  242.         if (gmtdefs.verbose) fprintf(stderr,"%d\t%d\t%.8lg\t%d\n", n_model, rank, c_chisq, 1);
  243.         if (robust) {
  244.             do {
  245.                 recompute_weights(data, n_data, work, &scale);
  246.                 move_model_a_to_b(c_model, w_model, n_model, &c_chisq, &w_chisq);
  247.                 load_gtg_and_gtd(data, n_data, gtg, gtd, workb, n_model, np);
  248.                 solve_system(gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, c_no, &rank);
  249.                 calc_m_and_r(data, n_data, c_model, n_model, workb);
  250.                 c_chisq = get_chisq(data, n_data, n_model);
  251.                 significant = sig_f(c_chisq, n_data-n_model, w_chisq, n_data-n_model, confid, &prob);
  252.                 if (gmtdefs.verbose) fprintf(stderr,"%d\t%d\t%.8lg\t%.3lg\n", n_model, rank, c_chisq, prob);
  253.             } while (significant);
  254.             /* Go back to previous model only if w_chisq < c_chisq  */
  255.             if (w_chisq < c_chisq) {
  256.                 move_model_a_to_b(w_model, c_model, n_model, &w_chisq, &c_chisq);
  257.                 calc_m_and_r(data, n_data, c_model, n_model, workb);
  258.                 if (weighted_output && n_model == n_model_max) recompute_weights(data, n_data, work, &scale);
  259.             }
  260.         }
  261.         /* First [robust] model has been found  */
  262.  
  263.         significant = TRUE;
  264.         while(n_model < n_model_max && significant) {
  265.             move_model_a_to_b(c_model, o_model, n_model, &c_chisq, &o_chisq);
  266.             n_model++;
  267.  
  268.             /* Fit next model  */
  269.             load_gtg_and_gtd(data, n_data, gtg, gtd, workb, n_model, np);
  270.             solve_system(gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, c_no, &rank);
  271.             calc_m_and_r(data, n_data, c_model, n_model, workb);
  272.             c_chisq = get_chisq(data, n_data, n_model);
  273.             if (gmtdefs.verbose) fprintf(stderr,"%d\t%d\t%.8lg\t%.3lg\n", n_model, rank, c_chisq, 1.00);
  274.             if (robust) {
  275.                 do {
  276.                     recompute_weights(data, n_data, work, &scale);
  277.                     move_model_a_to_b(c_model, w_model, n_model, &c_chisq, &w_chisq);
  278.                     load_gtg_and_gtd(data, n_data, gtg, gtd, workb, n_model, np);
  279.                     solve_system(gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, c_no, &rank);
  280.                     calc_m_and_r(data, n_data, c_model, n_model, workb);
  281.                     c_chisq = get_chisq(data, n_data, n_model);
  282.                     significant = sig_f(c_chisq, n_data-n_model, w_chisq, n_data-n_model, confid, &prob);
  283.                     if (gmtdefs.verbose) fprintf(stderr,"%d\t%d\t%.8lg\t%.3lg\n", n_model, rank, c_chisq, prob);
  284.                 } while (significant);
  285.                 /* Go back to previous model only if w_chisq < c_chisq  */
  286.                 if (w_chisq < c_chisq) {
  287.                     move_model_a_to_b(w_model, c_model, n_model, &w_chisq, &c_chisq);
  288.                     calc_m_and_r(data, n_data, c_model, n_model, workb);
  289.                     if (weighted_output && n_model == n_model_max) recompute_weights(data, n_data, work, &scale);
  290.                 }
  291.             }
  292.             /* Next [robust] model has been found  */
  293.             significant = sig_f(c_chisq, n_data-n_model, o_chisq, n_data-n_model-1, confid, &prob);
  294.         }
  295.  
  296.         if (!(significant) ) {    /* Go back to previous [robust] model, stored in o_model  */
  297.             n_model--;
  298.             rank--;
  299.             move_model_a_to_b(o_model, c_model, n_model, &o_chisq, &c_chisq);
  300.             calc_m_and_r(data, n_data, c_model, n_model, workb);
  301.             if (robust && weighted_output) recompute_weights(data, n_data, work, &scale);
  302.         }
  303.     }
  304.     else {
  305.         n_model = n_model_max;
  306.         load_gtg_and_gtd(data, n_data, gtg, gtd, workb, n_model, np);
  307.         solve_system(gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, c_no, &rank);
  308.         calc_m_and_r(data, n_data, c_model, n_model, workb);
  309.         c_chisq = get_chisq(data, n_data, n_model);
  310.         if (gmtdefs.verbose) fprintf(stderr,"%d\t%d\t%.8lg\t%.3lg\n", n_model, rank, c_chisq, 1.00);
  311.         if (robust) {
  312.             do {
  313.                 recompute_weights(data, n_data, work, &scale);
  314.                 move_model_a_to_b(c_model, w_model, n_model, &c_chisq, &w_chisq);
  315.                 load_gtg_and_gtd(data, n_data, gtg, gtd, workb, n_model, np);
  316.                 solve_system(gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, c_no, &rank);
  317.                 calc_m_and_r(data, n_data, c_model, n_model, workb);
  318.                 c_chisq = get_chisq(data, n_data, n_model);
  319.                 significant = sig_f(c_chisq, n_data-n_model, w_chisq, n_data-n_model, confid, &prob);
  320.                 if (gmtdefs.verbose) fprintf(stderr,"%d\t%d\t%.8lg\t%.3lg\n", n_model, rank, c_chisq, prob);
  321.             } while (significant);
  322.             /* Go back to previous model only if w_chisq < c_chisq  */
  323.             if (w_chisq < c_chisq) {
  324.                 move_model_a_to_b(w_model, c_model, n_model, &w_chisq, &c_chisq);
  325.                 calc_m_and_r(data, n_data, c_model, n_model, workb);
  326.                 if (weighted_output && n_model == n_model_max) recompute_weights(data, n_data, work, &scale);
  327.             }
  328.         }
  329.     }
  330.  
  331.     if (gmtdefs.verbose) {
  332.         fprintf(stderr,"Final model stats:  N model parameters %d.  Rank %d.  Chi-Squared:  %.8lg\n", n_model, rank, c_chisq);
  333.         fprintf(stderr,"Model Coefficients:");
  334.         for (i = 0; i < n_model; i++) {
  335.             fprintf(stderr, format, c_model[i]);
  336.         }
  337.         fprintf(stderr,"\n");
  338.     }
  339.  
  340.     untransform_x(data, n_data, xmin, xmax, ymin, ymax);
  341.  
  342.     write_output(data, n_data, output_choice, n_outputs);
  343.  
  344.     free_the_memory(gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
  345.  
  346.     gmt_end (argc, argv);
  347. }
  348.  
  349. void    read_data(data, n_data, xmin, xmax, ymin, ymax, weighted_input, work, fp)
  350. struct    DATA    **data;
  351. int    *n_data, weighted_input;
  352. double    *xmin, *xmax, *ymin, *ymax, **work;
  353. FILE    *fp;
  354. {
  355.  
  356.     int    i = 0, ix, iy, n_alloc = GMT_CHUNK;
  357.     char    buffer[512];
  358.     double    xy[2];
  359.  
  360.         if (fp == NULL) fp = stdin;
  361.         (*data) = (struct DATA *) memory (CNULL, n_alloc, sizeof(struct DATA), "trend1d");
  362.     ix = (gmtdefs.xy_toggle) ? 1 : 0;    iy = 1 - ix;        /* Set up which columns have x and y */
  363.  
  364.     while ( (fgets (buffer, 512, fp) ) != NULL) {
  365.         if ((sscanf(buffer, "%lf %lf %lf %lf", &xy[ix], &xy[iy], &(*data)[i].z, &(*data)[i].w)) != (3+weighted_input)) {
  366.             fprintf(stderr,"trend2d:  Cannot read x y z [w] from line %d\n", i+1);
  367.             continue;
  368.         }
  369.         (*data)[i].x = xy[0];
  370.         (*data)[i].y = xy[1];
  371.         if (!(weighted_input)) (*data)[i].w = 1.0;
  372.  
  373.         if (i) {
  374.             if (*xmin > (*data)[i].x) *xmin = (*data)[i].x;
  375.             if (*xmax < (*data)[i].x) *xmax = (*data)[i].x;
  376.             if (*ymin > (*data)[i].y) *ymin = (*data)[i].y;
  377.             if (*ymax < (*data)[i].y) *ymax = (*data)[i].y;
  378.         }
  379.         else {
  380.             *xmin = (*data)[i].x;
  381.             *xmax = (*data)[i].x;
  382.             *ymin = (*data)[i].y;
  383.             *ymax = (*data)[i].y;
  384.         }
  385.                 i++;
  386.  
  387.                 if (i == n_alloc) {
  388.                         n_alloc += GMT_CHUNK;
  389.                         *data = (struct DATA *) memory ((char *)*data, n_alloc, sizeof(struct DATA), "trend1d");
  390.                 }
  391.         }
  392.         if (fp != stdin) fclose(fp);
  393.         *data = (struct DATA *) memory ((char *)*data, i, sizeof(struct DATA), "trend1d");
  394.         *work = (double *) memory (CNULL, i, sizeof(double), "trend1d");
  395.     *n_data = i;
  396. }
  397.  
  398. void    allocate_array_space(np, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model)
  399. int    np;
  400. double    **gtg, **v, **gtd, **lambda, **workb, **workz, **c_model, **o_model, **w_model;
  401. {
  402.         *gtg = (double *) memory (CNULL, np*np, sizeof(double), "trend1d");
  403.         *v = (double *) memory (CNULL, np*np, sizeof(double), "trend1d");
  404.         *gtd = (double *) memory (CNULL, np, sizeof(double), "trend1d");
  405.         *lambda = (double *) memory (CNULL, np, sizeof(double), "trend1d");
  406.         *workb = (double *) memory (CNULL, np, sizeof(double), "trend1d");
  407.         *workz = (double *) memory (CNULL, np, sizeof(double), "trend1d");
  408.         *c_model = (double *) memory (CNULL, np, sizeof(double), "trend1d");
  409.         *o_model = (double *) memory (CNULL, np, sizeof(double), "trend1d");
  410.         *w_model = (double *) memory (CNULL, np, sizeof(double), "trend1d");
  411. }
  412.     
  413. void    write_output(data, n_data, output_choice, n_outputs)
  414. struct    DATA    *data;
  415. int    n_data, n_outputs;
  416. char    *output_choice;
  417. {
  418.     int    i, j;
  419.     char format1[20], format2[20];
  420.  
  421.     sprintf (format1, "%s\t\0", gmtdefs.d_format);
  422.     sprintf (format2, "%s\n\0", gmtdefs.d_format);
  423.  
  424.     for (i = 0; i < n_data; i++) {
  425.         for (j = 0; j < n_outputs-1; j++) {
  426.             switch (output_choice[j]) {
  427.                 case 'x':
  428.                     printf(format1, data[i].x);
  429.                     break;
  430.                 case 'y':
  431.                     printf(format1, data[i].y);
  432.                     break;
  433.                 case 'z':
  434.                     printf(format1, data[i].z);
  435.                     break;
  436.                 case 'm':
  437.                     printf(format1, data[i].m);
  438.                     break;
  439.                 case 'r':
  440.                     printf(format1, data[i].r);
  441.                     break;
  442.                 case 'w':
  443.                     printf(format1, data[i].w);
  444.                     break;
  445.             }
  446.         }
  447.         switch (output_choice[j]) {
  448.             case 'x':
  449.                 printf(format2, data[i].x);
  450.                 break;
  451.             case 'y':
  452.                 printf(format2, data[i].y);
  453.                 break;
  454.             case 'z':
  455.                 printf(format2, data[i].z);
  456.                 break;
  457.             case 'm':
  458.                 printf(format2, data[i].m);
  459.                 break;
  460.             case 'r':
  461.                 printf(format2, data[i].r);
  462.                 break;
  463.             case 'w':
  464.                 printf(format2, data[i].w);
  465.                 break;
  466.         }
  467.     }
  468. }
  469.  
  470. void    free_the_memory(gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work)
  471. double    *gtg, *v, *gtd, *lambda, *workb, *workz, *c_model, *o_model, *w_model, *work;
  472. struct    DATA    *data;
  473. {
  474.     free ( (char *)work);
  475.     free ( (char *)data);
  476.     free ( (char *)w_model);
  477.     free ( (char *)o_model);
  478.     free ( (char *)c_model);
  479.     free ( (char *)workz);
  480.     free ( (char *)workb);
  481.     free ( (char *)lambda);
  482.     free ( (char *)gtd);
  483.     free ( (char *)v);
  484.     free ( (char *)gtg);
  485. }
  486.  
  487. void    transform_x(data, n_data, xmin, xmax, ymin, ymax)
  488. struct    DATA    *data;
  489. int    n_data;
  490. double    xmin, xmax, ymin, ymax;
  491. {
  492.     int    i;
  493.     double    offsetx, scalex;
  494.     double    offsety, scaley;
  495.  
  496.     offsetx = 0.5 * (xmin + xmax);    /* Mid Range  */
  497.     offsety = 0.5 * (ymin + ymax);
  498.     scalex = 2.0 / (xmax - xmin);    /* 1 / (1/2 Range)  */
  499.     scaley = 2.0 / (ymax - ymin);
  500.  
  501.     for (i = 0; i < n_data; i++) {
  502.         data[i].x = (data[i].x - offsetx) * scalex;
  503.         data[i].y = (data[i].y - offsety) * scaley;
  504.     }
  505. }
  506.  
  507. void    untransform_x(data, n_data, xmin, xmax, ymin, ymax)
  508. struct    DATA    *data;
  509. int    n_data;
  510. double    xmin, xmax, ymin, ymax;
  511. {
  512.     int    i;
  513.     double    offsetx, scalex;
  514.     double    offsety, scaley;
  515.  
  516.     offsetx = 0.5 * (xmin + xmax);    /* Mid Range  */
  517.     offsety = 0.5 * (ymin + ymax);
  518.     scalex = 0.5 * (xmax - xmin);    /* 1/2 Range  */
  519.     scaley = 0.5 * (ymax - ymin);
  520.  
  521.     for (i = 0; i < n_data; i++) {
  522.         data[i].x = (data[i].x * scalex) + offsetx;
  523.         data[i].y = (data[i].y * scaley) + offsety;
  524.     }
  525. }
  526.  
  527. double    get_chisq(data, n_data, n_model)
  528. struct    DATA    *data;
  529. int    n_data, n_model;
  530. {
  531.     int    i, nu;
  532.     double    chi = 0.0;
  533.  
  534.  
  535.     for (i = 0; i < n_data; i++) {    /* Weight is already squared  */
  536.         if (data[i].w == 1.0) {
  537.             chi += (data[i].r * data[i].r);
  538.         }
  539.         else {
  540.             chi += (data[i].r * data[i].r * data[i].w);
  541.         }
  542.     }
  543.     nu = n_data - n_model;
  544.     if (nu > 1) return(chi/nu);
  545.     return(chi);
  546. }
  547.  
  548. void    recompute_weights(data, n_data, work, scale)
  549. struct    DATA    *data;
  550. int    n_data;
  551. double    *work, *scale;
  552. {
  553.     int    i;
  554.     double    k, ksq, rr;
  555.  
  556.     /* First find median { fabs(data[].r) },
  557.         estimate scale from this,
  558.         and compute chisq based on this.  */ 
  559.  
  560.     for (i = 0; i < n_data; i++) {
  561.         work[i] = fabs(data[i].r);
  562.     }
  563.     qsort( (char *)work, n_data, sizeof(double), comp_double_asc);
  564.  
  565.     if (n_data%2) {
  566.         *scale = 1.4826 * work[n_data/2];
  567.     }
  568.     else {
  569.         *scale = 0.7413 * (work[n_data/2 - 1] + work[n_data/2]);
  570.     }
  571.  
  572.     k = 1.5 * (*scale);    /*  Huber[1964] weight; 95% efficient for Normal data  */
  573.     ksq = k * k;
  574.  
  575.     for (i = 0; i < n_data; i++) {
  576.         rr = fabs(data[i].r);
  577.         if (rr <= k) {
  578.             data[i].w = 1.0;
  579.         }
  580.         else {
  581.             data[i].w = (2*k/rr) - (ksq/(rr*rr) );    /* This is really w-squared  */
  582.         }
  583.     }
  584. }
  585.  
  586. void    load_g_row(x, y, n, gr)
  587. double    x, y;    /* Current data position, appropriately normalized.  */
  588. int    n;    /* Number of model parameters, and elements of gr[]  */
  589. double    gr[];    /* Elements of row of G matrix.  */
  590. {
  591.     /* Routine computes the elements gr[j] in the ith row of the
  592.         G matrix (Menke notation), where x,y is the ith datum's
  593.         location.  */
  594.  
  595.     int    j;
  596.  
  597.     j = 0;
  598.     while (j < n) {
  599.         switch (j) {
  600.             case 0:
  601.                 gr[j] = 1.0;
  602.                 break;
  603.             case 1:
  604.                 gr[j] = x;
  605.                 break;
  606.             case 2:
  607.                 gr[j] = y;
  608.                 break;
  609.             case 3:
  610.                 gr[j] = x*y;
  611.                 break;
  612.             case 4:
  613.                 gr[j] = 2 * x * gr[1] - gr[0];
  614.                 break;
  615.             case 5:
  616.                 gr[j] = 2 * y * gr[2] - gr[0];
  617.                 break;
  618.             case 6:
  619.                 gr[j] = 2 * x * gr[4] - gr[1];
  620.                 break;
  621.             case 7:
  622.                 gr[j] = gr[4] * gr[2];
  623.                 break;
  624.             case 8:
  625.                 gr[j] = gr[5] * gr[1];
  626.                 break;
  627.             case 9:
  628.                 gr[j] = 2 * y * gr[5] - gr[2];
  629.                 break;
  630.         }
  631.         j++;
  632.     }
  633. }
  634.  
  635. void    calc_m_and_r(data, n_data, model, n_model, grow)
  636. struct    DATA    *data;
  637. int    n_data, n_model;
  638. double    *model, *grow;
  639. {
  640.     /*    model[n_model] holds solved coefficients of m_type model.
  641.         grow[n_model] is a vector for a row of G matrix.  */
  642.  
  643.     int    i, j;
  644.     for (i = 0; i < n_data; i++) {
  645.         load_g_row(data[i].x, data[i].y, n_model, grow);
  646.         data[i].m = 0.0;
  647.         for (j = 0; j < n_model; j++) {
  648.             data[i].m += model[j]*grow[j];
  649.         }
  650.         data[i].r = data[i].z - data[i].m;
  651.     }
  652. }
  653.  
  654. void    move_model_a_to_b(model_a, model_b, n_model, chisq_a, chisq_b)
  655. double    *model_a, *model_b, *chisq_a, *chisq_b;
  656. int    n_model;
  657. {
  658.     int    i;
  659.     for(i = 0; i<  n_model; i++) {
  660.         model_b[i] = model_a[i];
  661.     }
  662.     *chisq_b = *chisq_a;
  663. }
  664.  
  665. void    load_gtg_and_gtd(data, n_data, gtg, gtd, grow, n_model, mp)
  666. struct    DATA    *data;
  667. double    *gtg, *gtd, *grow;
  668. int    n_data, n_model, mp;    /* mp is row dimension of gtg  */
  669. {
  670.  
  671.     int    i, j, k;
  672.     double    wz;
  673.  
  674.     /* First zero the contents for summing:  */
  675.  
  676.     for (j = 0; j < n_model; j++) {
  677.         for (k = 0; k < n_model; k++) {
  678.             gtg[j + k*mp] = 0.0;
  679.         }
  680.         gtd[j] = 0.0;
  681.     }
  682.  
  683.     /* Sum over all data  */
  684.     for (i = 0; i < n_data; i++) {
  685.         load_g_row(data[i].x, data[i].y, n_model, grow);
  686.         if (data[i].w != 1.0) {
  687.             wz = data[i].w * data[i].z;
  688.             for (j = 0; j < n_model; j++) {
  689.                 for (k = 0; k < n_model; k++) {
  690.                     gtg[j + k*mp] += (data[i].w * grow[j] * grow[k]);
  691.                 }
  692.                 gtd[j] += (wz * grow[j]);
  693.             }
  694.         }
  695.         else {
  696.             for (j = 0; j < n_model; j++) {
  697.                 for (k = 0; k < n_model; k++) {
  698.                     gtg[j + k*mp] += (grow[j] * grow[k]);
  699.                 }
  700.                 gtd[j] += (data[i].z * grow[j]);
  701.             }
  702.         }
  703.     }
  704. }
  705.  
  706. void    solve_system(gtg, gtd, model, n_model, mp, lambda, v, b, z, c_no, ir)
  707. int    n_model, mp, *ir;
  708. double    *gtg, *gtd, *model, *lambda, *v, *b, *z, c_no;
  709.  
  710. {
  711.  
  712.     int    i, j, k, rank = 0, n, m, nrots;
  713.     double    c_test, temp_inverse_ij;
  714.  
  715.     if (n_model == 1) {
  716.         model[0] = gtd[0] / gtg[0];
  717.         *ir = 1;
  718.     }
  719.     else {
  720.         n = n_model;
  721.         m = mp;
  722.         if(jacobi(gtg, &n, &m, lambda, v, b, z, &nrots)) {
  723.             fprintf(stderr,"trend2d:  Warning:  Matrix Solver Convergence Failure.\n");
  724.         }
  725.         c_test = fabs(lambda[0])/c_no;
  726.         while(rank < n_model && lambda[rank] > 0.0 && lambda[rank] > c_test) rank++;
  727.         for (i = 0; i < n_model; i++) {
  728.             model[i] = 0.0;
  729.             for (j = 0; j < n_model; j++) {
  730.                 temp_inverse_ij = 0.0;
  731.                 for (k = 0; k <  rank; k++) {
  732.                     temp_inverse_ij += (v[i + k*mp] * v[j + k*mp] / lambda[k]);
  733.                 }
  734.                 model[i] += (temp_inverse_ij * gtd[j]);
  735.             }
  736.         }
  737.         *ir = rank;
  738.     }
  739. }
  740.