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