1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Brief synopsis: Reads stdin or file of x y pairs, or weighted pairs as x,y w data. Fit
19 * a regression model y = f(x) + e, where e are error misfits and f(x) has
20 * some user-prescribed functional form. Presently available models are
21 * polynomials and Fourier series. The user may choose the number of terms
22 * in the model to fit, whether to seek iterative refinement robust w.r.t.
23 * outliers, and whether to seek automatic discovery of the significant
24 * number of model parameters.
25 *
26 * Author: W. H. F. Smith
27 * Date: 1 JAN 2010
28 * Version: 6 API
29 *
30 * In trend1d I chose to construct the polynomial model using Chebyshev
31 * Polynomials so that the user may easily compare the sizes of the
32 * coefficients (and compare with a Fourier series as well). Tn(x)
33 * is an n-degree polynomial with n zero-crossings in [-1,1] and n+1
34 * extrema, at which the value of Tn(x) is +/- 1. It is this property
35 * which makes it easy to compare the size of the coefficients.
36 *
37 * During model fitting the data x coordinate is normalized into the domain
38 * [-1, 1] for Chebyshev Polynomial fitting, or into the domain [-pi, pi]
39 * for Fourier series fitting. Before writing out the data the coordinate
40 * is rescaled to match the original input values.
41 *
42 * An n degree polynomial can be written with terms of the form a0 + a1*x
43 * + a2*x*x + ... But it can also be written using other polynomial
44 * basis functions, such as a0*P0 + a1*P1 + a2*P2..., the Legendre
45 * polynomials, and a0*T0 + a1*T1 + a2*T2..., the Chebyshev polynomials.
46 * (The domain of the x values has to be in [-1, 1] in order to use P or T.)
47 * It is well known that the ordinary polynomial basis 1, x, x*x, ... gives
48 * terribly ill- conditioned matrices. The Ps and Ts do much better.
49 * This is because the ordinary basis is far from orthogonal. The Ps
50 * are orthogonal on [-1,1] and the Ts are orthogonal on [-1,1] under a
51 * simple weight function.
52 * Because the Ps have ordinary orthogonality on [-1,1], I expected them
53 * to be the best basis for a regression model; best meaning that they
54 * would lead to the most balanced G'G (matrix of normal equations) with
55 * the smallest condition number and the most nearly diagonal model
56 * parameter covariance matrix ((G'G)inverse). It turns out, however, that
57 * the G'G obtained from the Ts is very similar and usually has a smaller
58 * condition number than the Ps G'G. Both of these are vastly superior to
59 * the usual polynomials 1, x, x*x. In a test with 1000 equally spaced
60 * data and 8 model parameters, the Chebyshev system had a condition # = 10.6,
61 * Legendre = 14.8, and traditional = 54722.7. For 1000 randomly spaced data
62 * and 8 model parameters, the results were C = 13.1, L = 15.6, and P = 54916.6.
63 * As the number of model parameters approaches the number of data, the
64 * situation still holds, although all matrices get ill-conditioned; for 8
65 * random data and 8 model parameters, C = 1.8e+05, L = 2.6e+05, P = 1.1e+08.
66 * I expected the Legendre polynomials to have a covariance matrix more nearly
67 * diagonal than that of the Chebyshev polynomials, but on this criterion also
68 * the Chebyshev turned out to do better. Only as ndata -> n_model_parameters
69 * does the Legendre covariance matrix do better than the Chebyshev. So for
70 * all these reasons I use Chebyshev polynomials.
71 *
72 * Update: Aug-8, 2015 P. Wessel: Added ability to solve for a mixed model with
73 * both polynomial and Fourier parts. Also, ability to select just parts of a
74 * polynomial model (i.e., not include every term from 0 to n), but this
75 * necessitates working with powers of x and not Chebyshev, so we check and use
76 * the appropriate method.
77 */
78
79 #include "gmt_dev.h"
80
81 #define THIS_MODULE_CLASSIC_NAME "trend1d"
82 #define THIS_MODULE_MODERN_NAME "trend1d"
83 #define THIS_MODULE_LIB "core"
84 #define THIS_MODULE_PURPOSE "Fit [weighted] [robust] polynomial/Fourier model for y = f(x) to xy[w] data"
85 #define THIS_MODULE_KEYS "<D{,>D}"
86 #define THIS_MODULE_NEEDS ""
87 #define THIS_MODULE_OPTIONS "-:>Vbdefhiqsw" GMT_OPT("H")
88
89 #define TREND1D_N_OUTPUT_CHOICES 5
90
91 enum trend1d_enums {
92 TREND1D_NO_MODEL = 0,
93 TREND1D_POL_MODEL = 1,
94 TREND1D_POL_MODEL_NORM = 2,
95 TREND1D_CHEB_MODEL_NORM = 3
96 };
97
98 struct TREND1D_CTRL {
99 unsigned int n_outputs;
100 bool weighted_output;
101 unsigned int model_parameters; /* 0 = no output, 1 = polynomial output (users), 2 = polynomial output (normalized), 3 = Chebyshev (normalized) */
102 struct TREND1D_C { /* -C<condition_#> */
103 bool active;
104 double value;
105 } C;
106 struct TREND1D_F { /* -F<xymrw> */
107 bool active;
108 char col[TREND1D_N_OUTPUT_CHOICES]; /* Character codes for desired output in the right order */
109 } F;
110 struct TREND1D_I { /* -I[<confidence>] */
111 bool active;
112 double value;
113 } I;
114 struct TREND1D_N { /* -N[p|P|f|F|c|C|s|S|x|X]<list-of-terms>[,...][+l<length>][+o<origin>][+r] */
115 bool active;
116 struct GMT_MODEL M;
117 } N;
118 struct TREND1D_W { /* -W[+s] */
119 bool active;
120 unsigned int mode;
121 } W;
122 };
123
124 struct TREND1D_DATA {
125 double x; /* This is x and will be normalized */
126 double t; /* This will be in radians per x for Fourier terms */
127 double y;
128 double m;
129 double r;
130 double w;
131 };
132
trend1d_read_data(struct GMT_CTRL * GMT,struct TREND1D_DATA ** data,uint64_t * n_data,double * xmin,double * xmax,unsigned int weighted_input,double ** work)133 GMT_LOCAL int trend1d_read_data (struct GMT_CTRL *GMT, struct TREND1D_DATA **data, uint64_t *n_data, double *xmin, double *xmax, unsigned int weighted_input, double **work) {
134 uint64_t i;
135 size_t n_alloc = GMT_INITIAL_MEM_ROW_ALLOC;
136 double *in = NULL;
137 struct GMT_RECORD *In = NULL;
138
139 *data = gmt_M_memory (GMT, NULL, n_alloc, struct TREND1D_DATA);
140
141 i = 0;
142 do { /* Keep returning records until we reach EOF */
143 if ((In = GMT_Get_Record (GMT->parent, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
144 if (gmt_M_rec_is_error (GMT)) /* Bail if there are any read errors */
145 return (GMT_RUNTIME_ERROR);
146 else if (gmt_M_rec_is_any_header (GMT)) /* Skip all headers */
147 continue;
148 else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
149 break;
150 else
151 return (GMT_RUNTIME_ERROR); /* To shut up Coverity */
152 }
153
154 /* Data record to process */
155 in = In->data; /* Only need to process numerical part here */
156
157 (*data)[i].x = (*data)[i].t = in[GMT_X];
158 (*data)[i].y = in[GMT_Y];
159 if (weighted_input == 2) /* Got sigma, set weight = 1/s^2 */
160 (*data)[i].w = 1.0 / (in[GMT_Z] * in[GMT_Z]);
161 else if (weighted_input == 1) /* Got weight */
162 (*data)[i].w = in[GMT_Z];
163 else /* Default is unit weight */
164 (*data)[i].w = 1.0;
165
166 if (i) {
167 if (*xmin > (*data)[i].x) *xmin = (*data)[i].x;
168 if (*xmax < (*data)[i].x) *xmax = (*data)[i].x;
169 }
170 else {
171 *xmin = (*data)[i].x;
172 *xmax = (*data)[i].x;
173 }
174
175 if (++i == n_alloc) {
176 n_alloc <<= 1;
177 *data = gmt_M_memory (GMT, *data, n_alloc, struct TREND1D_DATA);
178 }
179 } while (true);
180
181 *data = gmt_M_memory (GMT, *data, i, struct TREND1D_DATA);
182 *work = gmt_M_memory (GMT, NULL, i, double);
183
184 *n_data = i;
185 return (0);
186 }
187
trend1d_allocate_the_memory(struct GMT_CTRL * GMT,unsigned int np,double ** gtg,double ** v,double ** gtd,double ** lambda,double ** workb,double ** workz,double ** c_model,double ** o_model,double ** w_model)188 GMT_LOCAL void trend1d_allocate_the_memory (struct GMT_CTRL *GMT, unsigned int np, double **gtg, double **v, double **gtd, double **lambda, double **workb, double **workz, double **c_model, double **o_model, double **w_model) {
189 *gtg = gmt_M_memory (GMT, NULL, np*np, double);
190 *v = gmt_M_memory (GMT, NULL, np*np, double);
191 *gtd = gmt_M_memory (GMT, NULL, np, double);
192 *lambda = gmt_M_memory (GMT, NULL, np, double);
193 *workb = gmt_M_memory (GMT, NULL, np, double);
194 *workz = gmt_M_memory (GMT, NULL, np, double);
195 *c_model = gmt_M_memory (GMT, NULL, np, double);
196 *o_model = gmt_M_memory (GMT, NULL, np, double);
197 *w_model = gmt_M_memory (GMT, NULL, np, double);
198 }
199
trend1d_write_output_trend(struct GMT_CTRL * GMT,struct TREND1D_DATA * data,uint64_t n_data,char * output_choice,unsigned int n_outputs)200 GMT_LOCAL void trend1d_write_output_trend (struct GMT_CTRL *GMT, struct TREND1D_DATA *data, uint64_t n_data, char *output_choice, unsigned int n_outputs) {
201 uint64_t i;
202 unsigned int j;
203 double out[5] = {0, 0, 0, 0, 0};
204 struct GMT_RECORD Out;
205
206 Out.data = out; Out.text = NULL;
207 for (i = 0; i < n_data; i++) {
208 for (j = 0; j < n_outputs; j++) {
209 switch (output_choice[j]) {
210 case 'x':
211 out[j] = data[i].x;
212 break;
213 case 'y':
214 out[j] = data[i].y;
215 break;
216 case 'm':
217 out[j] = data[i].m;
218 break;
219 case 'r':
220 out[j] = data[i].r;
221 break;
222 case 'w':
223 out[j] = data[i].w;
224 break;
225 }
226 }
227 GMT_Put_Record (GMT->parent, GMT_WRITE_DATA, &Out); /* Write this to output */
228 }
229 }
230
trend1d_free_the_memory(struct GMT_CTRL * GMT,double * gtg,double * v,double * gtd,double * lambda,double * workb,double * workz,double * c_model,double * o_model,double * w_model,struct TREND1D_DATA * data,double * work)231 GMT_LOCAL void trend1d_free_the_memory (struct GMT_CTRL *GMT, double *gtg, double *v, double *gtd, double *lambda, double *workb,
232 double *workz, double *c_model, double *o_model, double *w_model, struct TREND1D_DATA *data, double *work) {
233 gmt_M_free (GMT, work);
234 gmt_M_free (GMT, data);
235 gmt_M_free (GMT, w_model);
236 gmt_M_free (GMT, o_model);
237 gmt_M_free (GMT, c_model);
238 gmt_M_free (GMT, workz);
239 gmt_M_free (GMT, workb);
240 gmt_M_free (GMT, lambda);
241 gmt_M_free (GMT, gtd);
242 gmt_M_free (GMT, v);
243 gmt_M_free (GMT, gtg);
244 }
245
trend1d_transform_x(struct TREND1D_DATA * data,uint64_t n_data,struct GMT_MODEL * M,double xmin,double xmax)246 GMT_LOCAL void trend1d_transform_x (struct TREND1D_DATA *data, uint64_t n_data, struct GMT_MODEL *M, double xmin, double xmax) {
247 uint64_t i;
248 double offset, scale;
249
250 offset = (M->intercept) ? 0.5 * (xmin + xmax) : M->origin[GMT_X]; /* Mid Range or actual origin if no intercept */
251 scale = 2.0 / (xmax - xmin); /* 1 / (1/2 Range) */
252
253 /* Always normalize x for Chebyshev or polynomial fit */
254 if (M->type & 1) for (i = 0; i < n_data; i++) data[i].x = (data[i].x - offset) * scale; /* Now in -1/+1 range */
255 if (M->type & 2) { /* Have Fourier model component to deal with */
256 if (M->got_origin[GMT_X]) offset = M->origin[GMT_X]; /* Override the offset using given origin */
257 if (M->got_period[GMT_X]) scale = 2.0 / M->period[GMT_X]; /* Override the period using given period */
258 scale *= M_PI; /* Set Range to 1 period = 2 pi */
259 for (i = 0; i < n_data; i++) data[i].t = (data[i].t - offset) * scale; /* Now in units of period */
260 }
261 }
262
trend1d_untrend1d_transform_x(struct TREND1D_DATA * data,uint64_t n_data,struct GMT_MODEL * M,double xmin,double xmax)263 GMT_LOCAL void trend1d_untrend1d_transform_x (struct TREND1D_DATA *data, uint64_t n_data, struct GMT_MODEL *M, double xmin, double xmax) {
264 /* Undo transformation of x, if used */
265 uint64_t i;
266 double offset, scale;
267
268 if ((M->type & 1) == 0) return; /* Nothing to do */
269 offset = (M->intercept) ? 0.5 * (xmin + xmax) : M->origin[GMT_X]; /* Mid Range or actual origin if no intercept */
270 scale = 0.5 * (xmax - xmin); /* 1/2 Range */
271
272 for (i = 0; i < n_data; i++) data[i].x = (data[i].x * scale) + offset;
273 }
274
trend1d_get_chisq(struct TREND1D_DATA * data,uint64_t n_data,unsigned int n_model)275 GMT_LOCAL double trend1d_get_chisq (struct TREND1D_DATA *data, uint64_t n_data, unsigned int n_model) {
276 uint64_t i, nu;
277 double chi = 0.0;
278
279 for (i = 0; i < n_data; i++) { /* Weight is already squared */
280 if (data[i].w == 1.0)
281 chi += (data[i].r * data[i].r);
282 else
283 chi += (data[i].r * data[i].r * data[i].w);
284 }
285 nu = n_data - n_model;
286 if (nu > 1) return (chi/nu);
287 return (chi);
288 }
289
trend1d_recompute_weights(struct GMT_CTRL * GMT,struct TREND1D_DATA * data,uint64_t n_data,double * work,double * scale)290 GMT_LOCAL void trend1d_recompute_weights (struct GMT_CTRL *GMT, struct TREND1D_DATA *data, uint64_t n_data, double *work, double *scale) {
291 uint64_t i;
292 double k, ksq, rr;
293
294 /* First find median { fabs(data[].r) },
295 estimate scale from this,
296 and compute chisq based on this. */
297
298 for (i = 0; i < n_data; i++) work[i] = fabs(data[i].r);
299 gmt_sort_array (GMT, work, n_data, GMT_DOUBLE);
300
301 if (n_data%2)
302 *scale = MAD_NORMALIZE * work[n_data/2];
303 else
304 *scale = MAD_NORMALIZE * (work[n_data/2 - 1] + work[n_data/2]) / 2;
305
306 k = 1.5 * (*scale); /* Huber[1964] weight; 95% efficient for Normal data */
307 ksq = k * k;
308
309 for (i = 0; i < n_data; i++) {
310 rr = fabs(data[i].r);
311 data[i].w = (rr <= k) ? 1.0 : (2*k/rr) - (ksq/(rr*rr) ); /* This is really w-squared */
312 }
313 }
314
trend1d_chebyshev(double x,unsigned int n)315 GMT_LOCAL double trend1d_chebyshev (double x, unsigned int n) {
316 /* Return T_n(x) */
317 double cj, cj1, cj2;
318 unsigned int j;
319 if (n == 0) return 1.0;
320 if (n == 1) return x;
321 cj1 = 1.0; cj = x;
322 for (j = 2; j <= n; j++) {
323 cj2 = cj1;
324 cj1 = cj;
325 cj = 2.0 * x * cj1 - cj2;
326 }
327 return (cj);
328 }
329
trend1d_polynomial(double x,unsigned int n)330 GMT_LOCAL double trend1d_polynomial (double x, unsigned int n) {
331 /* Return x^n */
332 if (n == 0) return 1.0;
333 if (n == 1) return x;
334 return (pow (x, (double)n));
335 }
336
trend1d_load_g_row(double x,double t,int n,double * gr,struct GMT_MODEL * M)337 GMT_LOCAL void trend1d_load_g_row (double x, double t, int n, double *gr, struct GMT_MODEL *M) {
338 /* x: Current data position, appropriately normalized. */
339 /* Number of model parameters, and elements of gr[] */
340 /* Elements of row of G matrix. */
341 /* M: structure with info for each model term */
342
343 /* Routine computes the elements gr[j] in the ith row of the
344 G matrix (Menke notation), where x is the ith datum's
345 abscissa. */
346
347 int j, k;
348
349 for (j = 0; j < n; j++) {
350 k = M->term[j].order[GMT_X];
351 switch (M->term[j].kind) {
352 case GMT_POLYNOMIAL:
353 gr[j] = trend1d_polynomial (x, k);
354 break;
355 case GMT_CHEBYSHEV:
356 gr[j] = trend1d_chebyshev (x, k);
357 break;
358 case GMT_COSINE:
359 gr[j] = cos (k * t);
360 break;
361 case GMT_SINE:
362 gr[j] = sin (k * t);
363 break;
364 }
365 }
366 }
367
trend1d_calc_m_and_r(struct TREND1D_DATA * data,uint64_t n_data,double * model,unsigned int n_model,struct GMT_MODEL * M,double * grow)368 GMT_LOCAL void trend1d_calc_m_and_r (struct TREND1D_DATA *data, uint64_t n_data, double *model, unsigned int n_model, struct GMT_MODEL *M, double *grow) {
369 /* model[n_model] holds solved coefficients of m_type model.
370 grow[n_model] is a vector for a row of G matrix. */
371
372 uint64_t i;
373 unsigned int j;
374 for (i = 0; i < n_data; i++) {
375 trend1d_load_g_row (data[i].x, data[i].t, n_model, grow, M);
376 data[i].m = 0.0;
377 for (j = 0; j < n_model; j++) data[i].m += model[j]*grow[j];
378 data[i].r = data[i].y - data[i].m;
379 }
380 }
381
trend1d_move_model_a_to_b(double * model_a,double * model_b,unsigned int n_model,double * chisq_a,double * chisq_b)382 GMT_LOCAL void trend1d_move_model_a_to_b (double *model_a, double *model_b, unsigned int n_model, double *chisq_a, double *chisq_b) {
383 gmt_M_memcpy (model_b, model_a, n_model, double);
384 *chisq_b = *chisq_a;
385 }
386
trend1d_load_gtg_and_gtd(struct TREND1D_DATA * data,uint64_t n_data,double * gtg,double * gtd,double * grow,unsigned int n_model,unsigned int mp,struct GMT_MODEL * M)387 GMT_LOCAL void trend1d_load_gtg_and_gtd (struct TREND1D_DATA *data, uint64_t n_data, double *gtg, double *gtd, double *grow, unsigned int n_model, unsigned int mp, struct GMT_MODEL *M) {
388
389 /* mp is row dimension of gtg */
390
391 uint64_t i;
392 unsigned int j, k;
393 double wy;
394
395 /* First zero the contents for summing */
396
397 for (j = 0; j < n_model; j++) {
398 for (k = 0; k < n_model; k++) gtg[j + k*mp] = 0.0;
399 gtd[j] = 0.0;
400 }
401
402 /* Sum over all data */
403 for (i = 0; i < n_data; i++) {
404 trend1d_load_g_row (data[i].x, data[i].t, n_model, grow, M);
405 if (data[i].w != 1.0) {
406 wy = data[i].w * data[i].y;
407 for (j = 0; j < n_model; j++) {
408 for (k = 0; k < n_model; k++) gtg[j + k*mp] += (data[i].w * grow[j] * grow[k]);
409 gtd[j] += (wy * grow[j]);
410 }
411 }
412 else {
413 for (j = 0; j < n_model; j++) {
414 for (k = 0; k < n_model; k++) gtg[j + k*mp] += (grow[j] * grow[k]);
415 gtd[j] += (data[i].y * grow[j]);
416 }
417 }
418 }
419 }
420
trend1d_solve_system(struct GMT_CTRL * GMT,double * gtg,double * gtd,double * model,unsigned int n_model,unsigned int mp,double * lambda,double * v,double * b,double * z,double c_no,unsigned int * ir)421 GMT_LOCAL void trend1d_solve_system (struct GMT_CTRL *GMT, double *gtg, double *gtd, double *model, unsigned int n_model, unsigned int mp, double *lambda, double *v, double *b, double *z, double c_no, unsigned int *ir) {
422 unsigned int i, j, k, rank = 0, nrots;
423 double c_test, temp_inverse_ij;
424
425 if (n_model == 1) {
426 model[0] = gtd[0] / gtg[0];
427 *ir = 1;
428 }
429 else {
430 if (gmt_jacobi (GMT, gtg, n_model, mp, lambda, v, b, z, &nrots)) {
431 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Matrix Solver Convergence Failure.\n");
432 }
433 c_test = fabs (lambda[0]) / c_no;
434 while (rank < n_model && lambda[rank] > 0.0 && lambda[rank] > c_test) rank++;
435 for (i = 0; i < n_model; i++) {
436 model[i] = 0.0;
437 for (j = 0; j < n_model; j++) {
438 temp_inverse_ij = 0.0;
439 for (k = 0; k < rank; k++) {
440 temp_inverse_ij += (v[i + k*mp] * v[j + k*mp] / lambda[k]);
441 }
442 model[i] += (temp_inverse_ij * gtd[j]);
443 }
444 }
445 *ir = rank;
446 }
447 }
448
trend1d_unscale_polynomial(struct GMT_CTRL * GMT,double c[],unsigned int n,double a,double b)449 GMT_LOCAL void trend1d_unscale_polynomial (struct GMT_CTRL *GMT, double c[], unsigned int n, double a, double b) {
450 /* n are the first n terms that are polynomial - there may be Fourier terms beyond this set */
451 unsigned int j, k;
452 double cnst, fac;
453 gmt_M_unused(GMT);
454
455 cnst = fac = 2.0 / (b - a);
456 for (j = 1; j < n; j++) {
457 c[j] *= fac;
458 fac *= cnst;
459 }
460 if (n < 2) return; /* To avoid n-2 becoming huge since unsigned */
461 cnst = 0.5 * (a + b);
462 for (j = 0; j <= n - 2; j++) {
463 for (k = n - 1; k > j; k--) c[k-1] -= cnst * c[k];
464 }
465 }
466
trend1d_cheb_to_pol(struct GMT_CTRL * GMT,double c[],unsigned int un,double a,double b,unsigned int denorm)467 GMT_LOCAL void trend1d_cheb_to_pol (struct GMT_CTRL *GMT, double c[], unsigned int un, double a, double b, unsigned int denorm) {
468 /* Convert from Chebyshev coefficients used on a t = [-1,+1] interval
469 * to polynomial coefficients on the original x = [a b] interval.
470 * Modified from Numerical Miracles, ...eh Recipes */
471
472 int j, k, n = (int)un;
473 double sv, *d, *dd;
474
475 d = gmt_M_memory (GMT, NULL, n, double);
476 dd = gmt_M_memory (GMT, NULL, n, double);
477
478 /* First we generate coefficients for a polynomial in t */
479
480 d[0] = c[n-1];
481 for (j = n - 2; j >= 1; j--) {
482 for (k = n - j; k >= 1; k--) {
483 sv = d[k];
484 d[k] = 2.0 * d[k-1] - dd[k];
485 dd[k] = sv;
486 }
487 sv = d[0];
488 d[0] = -dd[0] + c[j];
489 dd[0] = sv;
490 }
491 for (j = n - 1; j >= 1; j--) d[j] = d[j-1] - dd[j];
492 /* d[0] = -dd[0] + 0.5 * c[0]; */ /* This is what Num. Rec. says, but we do not do the approx with 0.5 * c[0] */
493 d[0] = -dd[0] + c[0];
494
495 /* Next step is to undo the scaling so we can use coefficients with the user's x */
496
497 if (denorm)
498 trend1d_unscale_polynomial (GMT, d, un, a, b);
499
500 /* Return the new coefficients via c */
501
502 gmt_M_memcpy (c, d, un, double);
503
504 gmt_M_free (GMT, d);
505 gmt_M_free (GMT, dd);
506 }
507
New_Ctrl(struct GMT_CTRL * GMT)508 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
509 struct TREND1D_CTRL *C = NULL;
510
511 C = gmt_M_memory (GMT, NULL, 1, struct TREND1D_CTRL);
512
513 /* Initialize values whose defaults are not 0/false/NULL */
514 C->C.value = 1.0e06; /* Condition number for matrix solution */
515 C->I.value = 0.51; /* Confidence interval for significance test */
516 return (C);
517 }
518
Free_Ctrl(struct GMT_CTRL * GMT,struct TREND1D_CTRL * C)519 static void Free_Ctrl (struct GMT_CTRL *GMT, struct TREND1D_CTRL *C) { /* Deallocate control structure */
520 if (!C) return;
521 gmt_M_free (GMT, C);
522 }
523
usage(struct GMTAPI_CTRL * API,int level)524 static int usage (struct GMTAPI_CTRL *API, int level) {
525 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
526 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
527 GMT_Usage (API, 0, "usage: %s [<table>] -F<xymrw>|p|P|c -N[p|P|f|F|c|C|s|S|x|X]<list-of-terms>[,...][+l<length>][+o<origin>][+r] "
528 "[-C<condition_#>] [-I[<confidence>]] [%s] [-W[+s|w]] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
529 name, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, GMT_i_OPT, GMT_q_OPT, GMT_s_OPT, GMT_w_OPT, GMT_colon_OPT, GMT_PAR_OPT);
530
531 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
532
533 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
534 GMT_Option (API, "<");
535 GMT_Usage (API, -2, "Note: Input must provide (x,y[,w]) records; see -W for weights.");
536 GMT_Usage (API, 1, "\n-F<xymrw|p|P|c>");
537 GMT_Usage (API, -2, "Choose at least 1, up to 5, any order, of xymrw for output to standard output:");
538 GMT_Usage (API, 3, "x: The x-coordinate.");
539 GMT_Usage (API, 3, "y: The y-value.");
540 GMT_Usage (API, 3, "m: The model prediction.");
541 GMT_Usage (API, 3, "r: The residual (y-m).");
542 GMT_Usage (API, 3, "w: The weight (determined iteratively if robust fit used).");
543 GMT_Usage (API, -2, "Alternatively, request output of model coefficients instead:");
544 GMT_Usage (API, 3, "p: Polynomial coefficients.");
545 GMT_Usage (API, 3, "P: Normalized polynomial coefficients.");
546 GMT_Usage (API, 3, "c: Normalized Chebyshev coefficients.");
547 GMT_Usage (API, 1, "\n -N[p|P|f|F|c|C|s|S|x|X]<list-of-terms>[,...][+l<length>][+o<origin>][+r]");
548 GMT_Usage (API, -2, "Specify a polynomial, Fourier, or mixed model of any order; separate components by commas:");
549 GMT_Usage (API, 3, "p: Append <n> to add complete polynomial (including intercept) up to degree <n>.");
550 GMT_Usage (API, 3, "P: Append <n> (or xx..x) to add the component x^<n> only.");
551 GMT_Usage (API, 3, "f: Append <n> to add the Fourier series components 1-n.");
552 GMT_Usage (API, 3, "F: Append <n> to add just the <n>'th Fourier component.");
553 GMT_Usage (API, 3, "c: Append <n> to add the cosine series components 1-n.");
554 GMT_Usage (API, 3, "C: Append <n> to add just the <n>'th cosine component.");
555 GMT_Usage (API, 3, "s: Append <n> to add the sine series components 1-n.");
556 GMT_Usage (API, 3, "S: Append <n> to add just the <n>'th sine component.");
557 GMT_Usage (API, -2, "A few modifiers control further behavior:");
558 GMT_Usage (API, 3, "+l Append <period> to set fundamental period of x [range of x].");
559 GMT_Usage (API, 3, "+o Append <orig> to set origin of x [mid-point of x].");
560 GMT_Usage (API, 3, "+r Request robust model. E.g., robust quadratic = -Np2+r.");
561 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
562 GMT_Usage (API, 1, "\n-C<condition_#>");
563 GMT_Usage (API, -2, "Truncate eigenvalue spectrum so matrix has <condition_#> [Default = 1.0e06].");
564 GMT_Usage (API, 1, "\n-I[<confidence>]");
565 GMT_Usage (API, -2, "Iteratively Increase the number of model parameters, to a max of <n_model> so long as the "
566 "reduction in variance is significant at the <confidence> level [0.51].");
567 GMT_Option (API, "V");
568 GMT_Usage (API, 1, "\n-W[+s|w]");
569 GMT_Usage (API, -2, " Weighted input given, weights in 3rd column [Default is unweighted]. Select modifier:");
570 GMT_Usage (API, 3, "+s Read standard deviations and compute weights as 1/s^2.");
571 GMT_Usage (API, 3, "+w Read weights directly [Default].");
572 GMT_Option (API, "bi");
573 if (gmt_M_showusage (API)) GMT_Usage (API, -2, "Default is 2 (or 3 if -W is set) input columns.");
574 GMT_Option (API, "bo,d,e,f,h,i,q,s,w,:,.");
575
576 return (GMT_MODULE_USAGE);
577 }
578
parse(struct GMT_CTRL * GMT,struct TREND1D_CTRL * Ctrl,struct GMT_OPTION * options)579 static int parse (struct GMT_CTRL *GMT, struct TREND1D_CTRL *Ctrl, struct GMT_OPTION *options) {
580 /* This parses the options provided to trend1d and sets parameters in CTRL.
581 * Any GMT common options will override values set previously by other commands.
582 * It also replaces any file names specified as input or output with the data ID
583 * returned when registering these sources/destinations with the API.
584 */
585
586 unsigned int n_errors = 0, j;
587 struct GMT_OPTION *opt = NULL;
588 struct GMTAPI_CTRL *API = GMT->parent;
589
590 for (opt = options; opt; opt = opt->next) {
591 switch (opt->option) {
592
593 case '<': /* Skip input files */
594 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
595 break;
596
597 /* Processes program-specific parameters */
598
599 case 'C':
600 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
601 Ctrl->C.active = true;
602 Ctrl->C.value = atof (opt->arg);
603 break;
604 case 'F':
605 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
606 Ctrl->F.active = true;
607 for (j = 0; opt->arg[j]; j++) {
608 if (j < TREND1D_N_OUTPUT_CHOICES)
609 Ctrl->F.col[j] = opt->arg[j];
610 else {
611 GMT_Report (API, GMT_MSG_ERROR, "Option -F: Too many output columns selected: Choose from -Fxymrw|p\n");
612 n_errors++;
613 }
614 }
615 break;
616 case 'I':
617 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
618 Ctrl->I.active = true;
619 if (opt->arg[0]) Ctrl->I.value = atof (opt->arg);
620 break;
621 case 'N':
622 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
623 Ctrl->N.active = true;
624 if (gmt_parse_model (API->GMT, opt->option, opt->arg, 1U, &(Ctrl->N.M)) == -1) {
625 GMT_Report (API, GMT_MSG_ERROR, "Option -N: See usage for details.\n");
626 n_errors++;
627 }
628 break;
629 case 'W':
630 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
631 Ctrl->W.active = true;
632 if (gmt_validate_modifiers (GMT, opt->arg, 'W', "sw", GMT_MSG_ERROR)) n_errors++;
633 Ctrl->W.mode = (strstr (opt->arg, "+s")) ? 2 : 1;
634 break;
635
636 default: /* Report bad options */
637 n_errors += gmt_default_error (GMT, opt->option);
638 break;
639 }
640 }
641
642 n_errors += gmt_M_check_condition (GMT, Ctrl->C.value <= 1.0, "Option -C: Condition number must be larger than unity\n");
643 n_errors += gmt_M_check_condition (GMT, Ctrl->I.value < 0.0 || Ctrl->I.value > 1.0, "Option -C: Give 0 < confidence level < 1.0\n");
644 n_errors += gmt_M_check_condition (GMT, Ctrl->N.M.n_terms <= 0, "Option -N: A positive number of terms must be specified\n");
645 n_errors += gmt_check_binary_io (GMT, (Ctrl->W.active) ? 3 : 2);
646 for (j = Ctrl->n_outputs = 0; j < TREND1D_N_OUTPUT_CHOICES && Ctrl->F.col[j]; j++) {
647 if (!strchr ("xymrwpPc", Ctrl->F.col[j])) {
648 GMT_Report (API, GMT_MSG_ERROR, "Option -F: Unrecognized output choice %c\n", Ctrl->F.col[j]);
649 n_errors++;
650 }
651 else if (Ctrl->F.col[j] == 'w')
652 Ctrl->weighted_output = true;
653 else if (Ctrl->F.col[j] == 'p')
654 Ctrl->model_parameters = TREND1D_POL_MODEL;
655 else if (Ctrl->F.col[j] == 'P')
656 Ctrl->model_parameters = TREND1D_POL_MODEL_NORM;
657 else if (Ctrl->F.col[j] == 'c')
658 Ctrl->model_parameters = TREND1D_CHEB_MODEL_NORM;
659 Ctrl->n_outputs++;
660 }
661 n_errors += gmt_M_check_condition (GMT, Ctrl->n_outputs == 0, "Option -F: Must specify at least one output columns \n");
662 n_errors += gmt_M_check_condition (GMT, Ctrl->n_outputs > 1 && Ctrl->model_parameters,
663 "Option -F: When selecting model parameters, it must be the only output\n");
664
665 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
666 }
667
668 #define bailout(code) {gmt_M_free_options (mode); return (code);}
669 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
670
GMT_trend1d(void * V_API,int mode,void * args)671 EXTERN_MSC int GMT_trend1d (void *V_API, int mode, void *args) {
672 unsigned int i, n_model, rank, np;
673 int error = 0;
674 bool significant;
675
676 uint64_t n_data;
677
678 double *gtg = NULL, *v = NULL, *gtd = NULL, *lambda = NULL, *workb = NULL;
679 double *workz = NULL, *c_model = NULL, *o_model = NULL, *w_model = NULL, *work = NULL;
680 double xmin = 0.0, xmax = 0.0, c_chisq, o_chisq = 0.0, w_chisq, scale = 1.0, prob;
681
682 char format[GMT_BUFSIZ];
683
684 struct TREND1D_DATA *data = NULL;
685 struct TREND1D_CTRL *Ctrl = NULL;
686 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
687 struct GMT_OPTION *options = NULL;
688 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
689
690 /*----------------------- Standard module initialization and parsing ----------------------*/
691
692 if (API == NULL) return (GMT_NOT_A_SESSION);
693 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
694 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
695
696 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
697
698 /* Parse the command-line arguments */
699
700 if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
701 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
702 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
703 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
704
705 /*---------------------------- This is the trend1d main code ----------------------------*/
706
707 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
708 np = Ctrl->N.M.n_terms; /* Row dimension for matrices gtg and v */
709
710 if ((error = GMT_Set_Columns (API, GMT_IN, 2 + Ctrl->W.active, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
711 Return (error);
712 }
713 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_NONE, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
714 Return (API->error);
715 }
716 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
717 Return (API->error);
718 }
719 trend1d_allocate_the_memory (GMT, np, >g, &v, >d, &lambda, &workb, &workz, &c_model, &o_model, &w_model);
720 if ((error = trend1d_read_data (GMT, &data, &n_data, &xmin, &xmax, Ctrl->W.mode, &work)) != 0) {
721 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
722 Return (error);
723 }
724 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
725 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
726 Return (API->error);
727 }
728
729 if (xmin == xmax) {
730 GMT_Report (API, GMT_MSG_ERROR, "Min and Max value of input data are the same.\n");
731 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
732 Return (GMT_RUNTIME_ERROR);
733 }
734 if (n_data == 0) {
735 GMT_Report (API, GMT_MSG_ERROR, "Could not read any data.\n");
736 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
737 Return (GMT_RUNTIME_ERROR);
738 }
739 if (n_data < (uint64_t)Ctrl->N.M.n_terms) GMT_Report (API, GMT_MSG_ERROR, "Ill-posed problem; n_data < n_model_max.\n");
740
741 trend1d_transform_x (data, n_data, &(Ctrl->N.M), xmin, xmax); /* Set domain to [-1, 1] or [-pi, pi] */
742
743 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
744 sprintf (format,"Read %%" PRIu64 " data with X values from %s to %s\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
745 GMT_Report (API, GMT_MSG_INFORMATION, format, n_data, xmin, xmax);
746 GMT_Report (API, GMT_MSG_INFORMATION, "N_model%sRank%sChi_Squared%sSignificance\n", GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator);
747 }
748
749 sprintf (format, "%%d%s%%d%s%s%s%s\n", GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out);
750
751 if (Ctrl->I.active) {
752 n_model = 1;
753
754 /* Fit first model */
755 trend1d_load_gtg_and_gtd (data, n_data, gtg, gtd, workb, n_model, np, &(Ctrl->N.M));
756 trend1d_solve_system (GMT, gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, Ctrl->C.value, &rank);
757 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
758 c_chisq = trend1d_get_chisq (data, n_data, n_model);
759 GMT_Report (API, GMT_MSG_INFORMATION, format, n_model, rank, c_chisq, 1.0);
760 if (Ctrl->N.M.robust) {
761 do {
762 trend1d_recompute_weights (GMT, data, n_data, work, &scale);
763 trend1d_move_model_a_to_b (c_model, w_model, n_model, &c_chisq, &w_chisq);
764 trend1d_load_gtg_and_gtd (data, n_data, gtg, gtd, workb, n_model, np, &(Ctrl->N.M));
765 trend1d_solve_system (GMT, gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, Ctrl->C.value, &rank);
766 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
767 c_chisq = trend1d_get_chisq (data, n_data, n_model);
768 significant = gmt_sig_f (GMT, c_chisq, n_data-n_model, w_chisq, n_data-n_model, Ctrl->I.value, &prob);
769 GMT_Report (API, GMT_MSG_INFORMATION, format, n_model, rank, c_chisq, prob);
770 } while (significant);
771 /* Go back to previous model only if w_chisq < c_chisq */
772 if (w_chisq < c_chisq) {
773 trend1d_move_model_a_to_b (w_model, c_model, n_model, &w_chisq, &c_chisq);
774 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
775 if (Ctrl->weighted_output && n_model == Ctrl->N.M.n_terms) trend1d_recompute_weights (GMT, data, n_data, work, &scale);
776 }
777 }
778 /* First [robust] model has been found */
779
780 significant = true;
781 while (n_model < Ctrl->N.M.n_terms && significant) {
782 trend1d_move_model_a_to_b (c_model, o_model, n_model, &c_chisq, &o_chisq);
783 n_model++;
784
785 /* Fit next model */
786 trend1d_load_gtg_and_gtd (data, n_data, gtg, gtd, workb, n_model, np, &(Ctrl->N.M));
787 trend1d_solve_system (GMT, gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, Ctrl->C.value, &rank);
788 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
789 c_chisq = trend1d_get_chisq (data, n_data, n_model);
790 GMT_Report (API, GMT_MSG_INFORMATION, format, n_model, rank, c_chisq, 1.00);
791 if (Ctrl->N.M.robust) {
792 do {
793 trend1d_recompute_weights (GMT, data, n_data, work, &scale);
794 trend1d_move_model_a_to_b (c_model, w_model, n_model, &c_chisq, &w_chisq);
795 trend1d_load_gtg_and_gtd (data, n_data, gtg, gtd, workb, n_model, np, &(Ctrl->N.M));
796 trend1d_solve_system (GMT, gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, Ctrl->C.value, &rank);
797 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
798 c_chisq = trend1d_get_chisq (data, n_data, n_model);
799 significant = gmt_sig_f (GMT, c_chisq, n_data-n_model, w_chisq, n_data-n_model, Ctrl->I.value, &prob);
800 GMT_Report (API, GMT_MSG_INFORMATION, format, n_model, rank, c_chisq, prob);
801 } while (significant);
802 /* Go back to previous model only if w_chisq < c_chisq */
803 if (w_chisq < c_chisq) {
804 trend1d_move_model_a_to_b (w_model, c_model, n_model, &w_chisq, &c_chisq);
805 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
806 if (Ctrl->weighted_output && n_model == Ctrl->N.M.n_terms) trend1d_recompute_weights (GMT, data, n_data, work, &scale);
807 }
808 }
809 /* Next [robust] model has been found */
810 significant = gmt_sig_f (GMT, c_chisq, n_data-n_model, o_chisq, n_data-n_model-1, Ctrl->I.value, &prob);
811 }
812
813 if (!(significant) ) { /* Go back to previous [robust] model, stored in o_model */
814 n_model--;
815 rank--;
816 trend1d_move_model_a_to_b (o_model, c_model, n_model, &o_chisq, &c_chisq);
817 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
818 if (Ctrl->N.M.robust && Ctrl->weighted_output) trend1d_recompute_weights (GMT, data, n_data, work, &scale);
819 }
820 }
821 else {
822 n_model = Ctrl->N.M.n_terms;
823 trend1d_load_gtg_and_gtd (data, n_data, gtg, gtd, workb, n_model, np, &(Ctrl->N.M));
824 trend1d_solve_system (GMT, gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, Ctrl->C.value, &rank);
825 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
826 c_chisq = trend1d_get_chisq (data, n_data, n_model);
827 GMT_Report (API, GMT_MSG_INFORMATION, format, n_model, rank, c_chisq, 1.00);
828 if (Ctrl->N.M.robust) {
829 do {
830 trend1d_recompute_weights (GMT, data, n_data, work, &scale);
831 trend1d_move_model_a_to_b (c_model, w_model, n_model, &c_chisq, &w_chisq);
832 trend1d_load_gtg_and_gtd (data, n_data, gtg, gtd, workb, n_model, np, &(Ctrl->N.M));
833 trend1d_solve_system (GMT, gtg, gtd, c_model, n_model, np, lambda, v, workb, workz, Ctrl->C.value, &rank);
834 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
835 c_chisq = trend1d_get_chisq (data, n_data, n_model);
836 significant = gmt_sig_f (GMT, c_chisq, n_data-n_model, w_chisq, n_data-n_model, Ctrl->I.value, &prob);
837 GMT_Report (API, GMT_MSG_INFORMATION, format, n_model, rank, c_chisq, prob);
838 } while (significant);
839 /* Go back to previous model only if w_chisq < c_chisq */
840 if (w_chisq < c_chisq) {
841 trend1d_move_model_a_to_b (w_model, c_model, n_model, &w_chisq, &c_chisq);
842 trend1d_calc_m_and_r (data, n_data, c_model, n_model, &(Ctrl->N.M), workb);
843 if (Ctrl->weighted_output && n_model == Ctrl->N.M.n_terms) trend1d_recompute_weights (GMT, data, n_data, work, &scale);
844 }
845 }
846 }
847
848 /* Before output, convert back to polynomial coefficients, if desired */
849
850 if (Ctrl->model_parameters && Ctrl->N.M.n_poly) {
851 if (Ctrl->N.M.chebyshev) { /* Solved using Chebyshev, perhaps convert to polynomial coefficients */
852 if (Ctrl->model_parameters != TREND1D_CHEB_MODEL_NORM) {
853 char *kind[2] = {"user-domain", "normalized"};
854 GMT_Report (API, GMT_MSG_INFORMATION, "Convert from normalized Chebyshev to %s polynomial coefficients\n", kind[Ctrl->model_parameters-1]);
855 trend1d_cheb_to_pol (GMT, c_model, Ctrl->N.M.n_poly, xmin, xmax, Ctrl->model_parameters == TREND1D_POL_MODEL);
856 }
857 else
858 GMT_Report (API, GMT_MSG_INFORMATION, "Report normalized Chebyshev coefficients\n");
859 }
860 else if (Ctrl->model_parameters != TREND1D_POL_MODEL_NORM) {
861 GMT_Report (API, GMT_MSG_INFORMATION, "Convert from normalized polynomial to user-domain polynomial coefficients\n");
862 trend1d_unscale_polynomial (GMT, c_model, Ctrl->N.M.n_poly, xmin, xmax);
863 }
864 }
865
866 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
867 sprintf (format, "Final model stats: N model parameters %%d. Rank %%d. Chi-Squared: %s\n", GMT->current.setting.format_float_out);
868 GMT_Report (API, GMT_MSG_INFORMATION, format, n_model, rank, c_chisq);
869 if (!Ctrl->model_parameters) { /* Only give verbose feedback on coefficients if not requested as output */
870 if (Ctrl->N.M.type & 1) { /* Has polynomial component */
871 if (Ctrl->N.M.chebyshev)
872 trend1d_cheb_to_pol (GMT, c_model, Ctrl->N.M.n_poly, xmin, xmax, 1);
873 sprintf (format, "Model Coefficients (Polynomial");
874 }
875 if (Ctrl->N.M.type & 2) /* Has Fourier components */
876 strcat (format, " and Fourier");
877 strcat (format, "): ");
878 GMT_Report (API, GMT_MSG_INFORMATION, format);
879 sprintf (format, "%s%s", GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out);
880 for (i = 0; i < n_model; i++) GMT_Message (API, GMT_TIME_NONE, format, c_model[i]);
881 GMT_Message (API, GMT_TIME_NONE, "\n");
882 }
883 }
884
885 trend1d_untrend1d_transform_x (data, n_data, &(Ctrl->N.M), xmin, xmax);
886
887 i = (Ctrl->model_parameters) ? n_model : Ctrl->n_outputs;
888 if ((error = GMT_Set_Columns (API, GMT_OUT, i, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
889 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
890 Return (error);
891 }
892 if (GMT_Init_IO (GMT->parent, GMT_IS_DATASET, GMT_IS_NONE, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
893 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
894 Return (API->error);
895 }
896 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
897 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
898 Return (API->error);
899 }
900 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
901 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
902 Return (API->error);
903 }
904
905 if (!Ctrl->model_parameters) /* Write any or all of the 'xymrw' */
906 trend1d_write_output_trend (GMT, data, n_data, Ctrl->F.col, Ctrl->n_outputs);
907 else { /* Write only the model parameters */
908 struct GMT_RECORD Rec;
909 Rec.data = c_model; Rec.text = NULL;
910 GMT_Put_Record (API, GMT_WRITE_DATA, &Rec);
911 }
912
913 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
914 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
915 Return (API->error);
916 }
917
918 trend1d_free_the_memory (GMT, gtg, v, gtd, lambda, workb, workz, c_model, o_model, w_model, data, work);
919
920 Return (GMT_NOERROR);
921 }
922